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C  cylinder  diameter  (m) 

Cj  the  control  point  on  the  /th  line  segment  (control  element),  approximating  the  air¬ 
foil  surface 

CD  drag  coefficient  (dimensionless) 

E  total  collision  efficiency  (%) 

Em  maximum  local  collision  efficiency  (%) 

g  acceleration  due  to  gravity  (m  s-2) 

K  Langmuir  inertia  parameter  (dimensionless) 

L  nondimensional  distance  along  the  surface  of  the  accretion,  starting  at  the  nose 
(dimensionless) 

N  number  of  line  segments  (control  elements)  approximating  the  airfoil  surface 
P  any  point  in  the  airstream 

P  air  pressure  (Pa) 

r  distance  between  a  point  on  a  control  element  and  any  point  in  the  airstream  (m) 

r(L)  local  radius  of  curvature  of  the  accretion  or  substrate  at  distance  L  from  the  nose 

(m) 

rd  droplet  radius  (m) 

R(L)  icing  flux  at  distance  L  from  the  nose  (kg  nr2  s’1) 

Re  Reynolds  number  of  the  droplet  (dimensionless) 

Sj  any  point  on  the  jth  control  element 
T  air  temperature  (K) 

f  time  (s) 

/A  total  accretion  time  for  a  layer  (s) 

u  x  component  of  airspeed  (m  s"') 

v  v  component  of  airspeed  (m  s’1) 

vx  ,v  component  of  droplet  impact  speed  (m  s  ') 

v  y  component  of  droplet  impact  speed  (m  s'1) 

V  ^  freest  ream  airspeed  (m  s'1) 

V,d  vector  air  velocity  (m  s’1) 

fd  vector  droplet  velocity  (m  s"') 

w  liquid  water  content  of  cloud  (kg  m"') 
v  A*-coordinate  (m) 

X  nondimensional  .v-coordinate  =  x/C  (dimensionless) 

Ad  nondimensional  droplet  position  vector  (dimensionless) 
v  .^-coordinate  (m) 

Y  nondimensional  v-coordinate  =  v/C  (dimensionless) 
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/30  maximum  local  collision  efficiency  (Vo) 
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gj  density  of  accreted  ice  (kg  m'3) 
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r  time  (s) 
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COMPUTER  MODELING  OF 
TIME-DEPENDENT  RIME  ICING 
IN  THE  ATMOSPHERE 


Edward  P.  Lozowski  and  Myron  M.  Oleskiw 


INTRODUCTION 

The  literature  on  the  subject  of  icing  is  very  extensive,  and  we  do  not  intend  to  review  it 
here.  Instead,  we  will  mention  simply  that  the  present  work  arose  chiefly  as  a  result  of  two 
earlier  investigations  into  icing,  one  at  the  U.S.  Army  Cold  Regions  Research  and  Engineer¬ 
ing  Laboratory  and  the  other  at  the  National  Research  Council  of  Canada  in  Ottawa.  The 
first  of  these  studies  was  reported  by  Ackley  and  Templeton  (1979),  while  the  second  was  de¬ 
scribed  by  Lozowski  et  al.  (1979).  Both  were  computer-simulated  models  of  ice  accretion  on 
a  cylinder.  The  first  included  time-dependent  effects  but  ignored  runback,  and  the  second  ig¬ 
nored  time  dependence  but  allowed  for  the  thermodynamics  of  runback. 

Although  cylinder  icing  models  are  of  intrinsic  importance  for  understanding  powerline  ic¬ 
ing,  for  example,  their  geometry  is  not  appropriate  for  the  study  of  airfoil  icing.  Airfoil  icing 
has  been  a  subject  of  renewed  interest  in  recent  years,  in  part  because  of  a  need  to  certify 
helicopters  and  general  aviation  aircraft  for  flight  in  IFR  (instrument  flight  rules)  icing  con¬ 
ditions.  The  limited  power  available  on  such  aircraft  and  the  new  materials  used  in  airfoil 
construction  demand  that  deicing  or  anti-icing  equipment  be  carefully  designed  for  maxi¬ 
mum  efficiency.  Although  the  design  of  such  equipment  requires  wind-tunnel  and  ultimately 
field  testing,  computer  simulation  models  are  considered  to  be  an  important  tool  in  the  de¬ 
sign  process  (Rosen  and  Potash  1981). 

The  objective  of  the  present  work  is  therefore  to  develop  and  test  a  computer  simulation 
model  for  airfoil  icing.  This  report  describes  a  model  that  permits  simulation  of  the  time-de¬ 
pendent  growth  of  ice  without  runback  on  an  arbitrary,  two-dimensional  airfoil.  In  develop¬ 
ing  the  model,  a  great  deal  of  effort  has  gone  into  carefully  specifying  the  assumptions  made 
and  into  testing  the  individual  components  of  the  model.  We  are  confident  that  within  the 
framework  of  assumptions  of  the  model,  the  icing  accretions  that  it  predicts  are  believable. 
Because  of  the  effort  required  to  develop  and  test  the  present  model,  it  has  not  been  possible 
in  the  time  available  to  make  the  model  completely  general.  Consequently,  we  have  not,  for 
example,  incorporated  accretion  thermodynamics  into  the  model  nor  taken  into  account 
rotation  effects,  such  as  could  be  found  on  helicopter  rotor  blades.  These  are  developments 
for  the  future.  Nevertheless,  the  model  as  it  stands  should  be  very  useful  for  estimating  the 
icing  rate  and  shape  on  airfoils  when  the  accretion  is  dry  (i.c.  no  runback)  and  when  rotation 
effects  (or  any  other  three-dimensional  effects)  may  be  ignored.  During  the  course  of  this 
work,  two  opportunties  arose  to  make  presentations  of  the  progress  to  date  to  audiences  of 
cloud  physicists  and  aerodynamicists.  These  presentations  are  summarized  in  Oleskiw'  and 
Lozowski  (1980)  and  Lozowski  and  Oleskiw'  (1981). 


METHODOLOGY 


The  modeling  of  airfoil  icing  may  be  separated  into  two  distinct  aspects.  The  first  is  the 
impingement  of  supercooled  water  droplets  on  the  airfoil  surface,  and  the  second  is  the  me¬ 
chanics  and  thermodynamics  of  the  resulting  accretion.  The  present  study  deals  exclusively 
with  the  first  aspect.  This  is  sufficient  for  the  investigation  of  the  dry  accretion  process,  in 
which  the  heat  transfer  is  great  enough  that  all  of  the  impinging  water  droplets  freeze  at  their 
point  of  impact.  This  restriction  is  analogous  to  that  made  by  Ackley  and  Templeton  (1979) 
in  their  model  of  icing  cylinders.  Cansdale  and  McNaughtan  (1977)  and  Lozowski  et  al. 
(1979)  also  considered  the  case  of  wet  accretion  on  cylinders,  in  which  some  impacting  water 
remains  unfrozen  and  is  blown  back  along  the  icing  surface.  The  extension  of  these  wet  ac¬ 
cretion  models  to  airfoils  requires  an  ability  to  calculate  the  heat  transfer  of  iced  airfoils.  We 
have  not  had  the  time  or  funds  to  do  this,  either  theoretically  or  by  experiment,  under  the 
present  contract. 

The  computer  algorithm  for  simulating  the  dry  accretion  process  may  be  broken  down  in¬ 
to  the  following  steps: 

1.  Determining  the  potential  flow  stream  function  field  around  an  arbitrary  two- 
dimensional  airfoil  in  crossflow. 

2.  Determining  the  incompressible  velocity  field  around  the  airfoil. 

3.  Calculating  droplet  trajectories  and  points  of  impact. 

4.  Determining  the  airfoil  collision  efficiency  as  a  function  of  surface  position  for 
specified  values  of  freestream  airspeed,  droplet  size,  and  airfoil  angle  of  attack. 

5.  Calculating  the  spatial  distribution  of  icing  during  a  short  time  interval,  under 
the  dry  accretion  assumption. 

6.  Determining  the  accretion  shape  and  mass. 

7.  Calculating  the  new  airfoil  shape  as  modified  by  the  ice  accretion. 

8.  Repeating  steps  one  to  seven  as  often  as  desired  to  obtain  the  growth  of  the  ac¬ 
cretion  as  a  function  of  time. 

In  the  detailed  descriptions  that  follow,  we  deal  with  each  of  these  steps  in  turn. 

Potential  flow  around  an  arbitrary  airfoil 

There  are  numerous  potential  flow  codes  available  that  permit  the  determination  of  the 
stream  function  field  around  an  arbitrary  two-dimensional  airfoil.  These  can  be  broadly 
classified  into  two  groups:  a)  conformal  transformation  techniques,  e.g.  Theodorson  and 
Garrick  (1932),  and  b)  surface  singularity  methods,  e.g.  Hess  and  Smith  (1967).  The  particu¬ 
lar  technique  chosen  to  address  the  icing  problem  should  have  the  following  characteristics. 
First,  it  should  be  particularly  efficient  in  terms  of  computer  time,  because  of  the  large 
number  of  air  velocity  calculations  required  to  determine  droplet  trajectories.  Secondly,  it 
should  be  capable  of  handling  the  changes  to  the  airfoil  profile  due  to  the  ice  accretion.  In 
this  latter  connection,  the  computer  code  must  be  capable  of  accepting  a  specification  of  the 
airfoil  in  terms  of  surface  coordinates  and,  moreover,  it  should  not  be  too  sensitive  to  small 
errors  in  the  specified  airfoil  coordinates. 

In  keeping  with  these  considerations,  we  chose  the  method  described  by  Kennedy  and 
Marsden  (1976).  This  is  one  of  the  so-called  “surface  singularity”  or  “panel”  methods.  It  is 
thought  to  be  the  simplest  available  and  provides  exceptional  accuracy  for  little  computing 
effort.  For  the  purpose  of  calculating  the  potential  flow,  the  airfoil  surface  is  approximated 
by  N  straight-line  segments  or  “panels,”  labeled  S-,  j  =  1,2,  .  .  .  ,  /V.  A  constant,  but  un¬ 
known,  vorticity  density  >(,Sj)  is  distributed  along  each  panel  or  control  element.  If  this  air¬ 
foil  model  is  immersed  in  a  uniform  stream  of  unit  velocity  (nondimensionalized),  at  an 
angle  of  attack  the  stream  function  at  any  point  external  to  the  airfoil  P(.v, y)  is,  according 
to  elementary  potential  How  theory,  given  by: 
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(1) 


=  y  cosa-x  sina  “  2“  53  /s.  lnr(P,5j)  r/Sj 

where  r(Pt5J)  is  the  distance  between  point  P  and  any  point  on  the  element  Sy 

To  solve  for  the  unknown  7(5j)  for  each  panel,  eq  1  is  applied  at  a  control  point,  C,,  i  = 

3,2 . V,  on  each  panel.  Imposing  the  boundary  condition,  that  the  stream  function  be  a 

constant  alu.ig  the  airfoil  surface  (i.e.  at  each  control  point),  and  the  Kutta  condition,  that 
the  surface  sti  camline  leave  the  trailing  edge  smoothly,  leads  to  a  set  of  linear  algebraic  equa¬ 
tions  for  the  unknown  vorticity  densities.  These  matrix  equations  are  solved  in  the  usual 
way,  and  eq  1  then  allows  the  determination  of  the  stream  function  anywhere  in  the  potential 
flow.  The  airspeed  at  any  point  can  then  be  determined  by  differentiation  of  eq  l. 

Other  investigators  (Bragg  and  Gregorek  1981)  have  adopted  the  conformal  transforma¬ 
tion  approach,  while  still  others  (McComber  and  Touzot  1981)  have  solved  Poisson’s  equa¬ 
tion  for  the  stream  function  using  finite  element  methods.  We  believe,  however,  that  the 
present  method  yields  greater  accuracy  and  spatial  resolution  for  a  similar  computing  effort. 

Incompressible  velocity  field 

The  air  velocity  components  may  be  calculated  at  any  desired  point  in  the  airstream  by  dif¬ 
ferentiating  the  stream  function;  that  is,  by  approximating  the  equations: 
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with  f  inite  differences.  This  is  done  with  a  space  increment,  Ax  or  Ay,  equal  to  the  diameter 
of  a  cloud  droplet .  Thus  it  is  possible  to  obtain  a  very  accurate  estimate  of  the  airspeed  at  the 
position  of  the  droplet’s  center  of  mass,  wherever  that  happens  to  be.  We  believe  that  this 
approach  is  more  accurate  than  that  used  by  some  other  investigators  (e.g.  Cansdale  1980, 
private  communication),  who  determine  the  airstream  velocity  field  initially  at  a  fixed  array 
of  grid  points.  When  the  air  velocity  at  the  droplet  position  is  desired,  an  interpolation  pro¬ 
cedure  among  the  grid  point  values  is  applied.  Our  approach  of  evaluating  the  air  velocity  as 
needed  at  precise  points  along  the  droplet  trajectory,  rather  than  by  interpolation,  is  made 
possible  by  the  economy  with  which  \fr  can  be  calculated  using  the  Kennedy-Marsden  ap¬ 
proach. 


Droplet  trajectory  equation 

Pearcey  and  Hill  (1956)  have  expressed  the  equation  of  motion  of  a  spherical  droplet  of 
fixed  mass  in  an  accelerated  air  flow  as: 
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where  A'd(jrd,.yd) 
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droplet  position  vector 

droplet  velocity  vector 

air  velocity  vector 

gravitational  acceleration 

steady-state  droplet  drag  coefficient 

air  density 

droplet  density 

droplet  radius 

kinematic  viscosity  of  the  air 
time. 


The  first  term  on  the  righthand  side  of  eq  4  is  the  net  buoyancy  of  the  droplet  in  air.  The 
second  term  is  the  steady  drag,  and  the  third  is  known  as  the  history  term  (because  of  the 
time  integral  over  the  entire  droplet  history).  The  first  two  terms  are  probably  in  need  of  no 
explanation,  although  the  gravitational  term  is  frequently  ignored  in  icing  calculations  (e.g. 
Langmuir  and  Blodgett  1946).  The  significance  of  the  history  term,  however,  may  not  be  so 
apparent.  It  is  essentially  a  correction  to  the  drag  term,  which  is  necessary  when  the  drag  co¬ 
efficient  used  in  the  second  term  is  the  steady-state  value,  appropriate  for  nonaccelerating 
droplets.  For  a  given  relative  velocity  between  the  droplet  and  the  airstream,  the  true  value  of 
CD  is  smaller  for  a  drop  that  is  accelerating  with  respect  to  the  flow  than  for  one  that  is  not 
accelerating  (i.e.  one  that  is  in  equilibrium  or  steady  state).  This  may  be  thought  of  as  a 
phase  lag  effect,  due  to  the  finite  rate  of  vorticity  diffusion,  which  requires  a  certain  time  for 
the  droplet  to  reach  equilibrium  with  the  airstream.  Because  of  the  large  droplet  acceleration 
that  occurs  in  certain  icing  situations,  we  felt  it  important  to  examine  the  effects  of  the  his¬ 
tory  term  on  the  calculation  of  the  droplet  trajectories.  Consequently,  comparisons  have 
been  made  between  results  calculated  without  the  history  term  (referred  to  as  the  steady-state 
drag  formulation)  and  those  calculated  with  the  history  term  included  (referred  to  as  the  non¬ 
steady-state  formulation).  It  should  be  noted  that  eq  4  also  incorporates  the  effects  of  the 
droplet’s  induced  mass  resulting  from  the  momentum  it  imparts  to  the  air  as  it  accelerates. 

The  formulation  used  to  determine  the  steady-state  drag  coefficient  as  a  function  of  drop¬ 
let  Reynolds  number  is  given  below: 


1.  Re  <0.01  CD  =  24/Red 

2.  0.01  <  Re  <5  CD=  24/Red  +  2.2 

3.  5  <  Re  <  5000  CD  =  0.2924(1 +9.06  Red0  5)2- 

The  droplet  Reynolds  number  is  defined  by: 

2rd63 


(5) 


Red  = 


where  and  f'a  are  respectively  the  droplet  and  air  velocity  vectors,  and  na  is  the  dynamic 
viscosity  of  the  air.  The  second  formulation  is  from  Sartor  and  Abbott  (1975),  while  the 
third  is  given  by  Abraham  (1970). 


Computational  procedure  for  trajectories 

Equations  3  and  4  in  component  form  yield  four  equations  that,  for  the  steady-state  for¬ 
mulation,  are  numerically  integrated  using  a  fourth-order  Runge-Kutta-Fehlberg  method 
(Lapidus  and  Seinfeld  1971,  Burden  et  ai.  1978).  This  procedure  permits  the  time  step  to  be 
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adjusted  continuously  for  optimum  speed  of  computation  given  a  specified  degree  of  accu¬ 
racy  required. 

When  the  history  term  is  incorporated  into  eq  4,  it  becomes  a  Voltera  integro-differential 
equation  of  the  second  kind.  The  method  of  solution  we  used  in  this  case  is  essentially  the 
same  as  that  used  for  the  steady-state  case,  with  the  additional  provision  that  the  history 
term  is  approximated  by  a  combined  numerical  and  analytical  technique.  With  this  scheme, 
the  integral  is  approximated  by  a  finite  sum  at  full  time  steps  of  the  Runge-Kutta-Fehlberg 
method  (for  example,  at  r  and  r+  At).  At  intermediate  time  steps,  however  (between  r  and 
t  +  At),  the  value  of  the  integral  is  approximated  by  the  extrapolation  of  a  Legendre  polyno¬ 
mial  fitted  to  the  previous  values  of  the  integral  at  full  time  steps. 


Determining  the  point  of  impact 

The  droplet  is  assumed  to  have  impinged  upon  the  airfoil  if  any  part  of  it  contacts  the  air¬ 
foil  surface.  Thus,  close  to  the  airfoil,  the  finite  size  of  the  droplet  is  taken  into  account.  This 
is  particularly  important  for  those  trajectories  just  within  the  envelope  of  colliding  trajecto¬ 
ries,  where  the  angle  of  incidence  from  the  normal  to  the  airfoil  surface  is  dose  to  90°. 

Calculation  of  collision  efficiencies 

To  determine  the  local  collision  efficiency,  /tf,  at  any  point  on  the  airfoil  surface,  use  is 
made  of  the  relation: 

,/#x  dY 

&(L)  =  -77  COSC* 


where  Y  is  the  ordinate  at  the  starting  point  of  a  particular  trajectory,  L  is  the  distance  along 
the  airfoil  surface  between  the  nose  and  the  impact  point  of  the  same  trajectory,  and  a  is  the 
angle  of  attack.  By  calculating  the  trajectories  for  between  10  and  20  droplets ,  Y  may  be 
plotted  as  a  function  of  Ly  and  the  derivative  taken  to  obtain  /i.  These  latter  operations  are  in 
fact  performed  numerically  using  a  cubic  spline  fit  to  the  point  values  on  the  graph  of  Y vs  L. 
The  resulting  local  collision  efficiencies  may  be  plotted  as  a  function  of  L  (e.g.  Fig.  2b)  or  of 
the  corresponding  abscissa  X  (e.g.  f  ig.  3a). 


Accreting  an  ice  layer 

In  the  model,  it  is  assumed  that  the  ice  growth  on  a  particular  small  segment  of  the  airfoil 
surface  is  oriented  perpendicular  to  the  surface.  According  to  Lozowski  et  al.  (1979),  the  ac¬ 
cretion  thickness  is  then  given  by  the  equation: 


h(L)  = 


2  R(L)tA 
Qf  U) 


(7) 


where 

R(L)  =  (8) 

is  the  icing  flux  with  1  ^  the  freestream  velocity  and  \v  the  liquid  water  content  of  the  air- 
stream.  tiSi  is  the  period  of  accretion,  q]  the  assumed  ice  density  (890  kg  nr'),  and  Mhc  radius 
of  curvature  of  the  airfoil  surface. 

In  the  results  presented  here,  we  assume  that  time  interval  is  sufficiently  small  that  the 
second  term  under  the  root  in  the  denominator  may  be  ignored. 

By  plotting  the  accretion  thickness  as  a  function  of  distance  along  the  airfoil  surface  from 
the  nose,  it  is  possible  to  determine  a  new  airfoil  surface  shape  after  it  has  iced  for  the  speci- 


fied  time  interval.  The  entire  procedure  can  now  be  repeated,  using  the  new  iced  airfoil  sur¬ 
face  to  determine  a  new  stream  function,  new  droplet  trajecu  >s,  and  ultimately  a  second 
accretion  layer.  By  continuing  in  this  manner,  it  is  possible  to  build  up  a  substantial  ice  accre¬ 
tion  on  the  airfoil. 


Determining  the  accuracy  of  the  flow  field 

The  accuracy  ot  the  Kennedy-Marsden  technique  was  tested  by  comparing  its  predicted 
stream  function  for  potential  flow  around  a  cylinder  with  the  known  analytic  solution.  Using 
50  control  elements,  the  error  in  0  is  at  most  0.1%  near  the  cylinder.  It  falls  to  below  0.01  % 
at  distances  from  the  surface  exceeding  about  four  cylinder  diameters. 

We  have  also  made  a  comparison  between  the  corresponding  velocity  fields.  In  one  such 
test  for  example,  using  a  cylinder  diameter  of  0.15  m,  an  air  pressure  of  78.5  kPa,  an  air  tem¬ 
perature  of  -  10°C,  and  a  freestream  velocity  of  1 14.3  m  s  ',  the  velocity  field  of  the  analytic 
solution  was  compared  with  that  provided  by  the  model  using  38  control  elements.  At  a 
distance  of  five  diameters  upstream,  the  air  velocities  differed  by  less  than  10  Very  close 
to  the  cylinder  they  were  as  high  as  1  to  2%.  However,  the  effect  of  these  airstream  velocity 
errors  on  the  computed  droplet  collision  efficiencies  was  found  to  be  much  less  than  1%. 


Determining  the  accuracy  of  the  trajectories 

To  establish  some  confidence  in  the  trajectories  themselves,  it  was  decided  to  make  a  com¬ 
parison  with  two  cases  considered  by  Langmuir  and  Blodgett  (1946).  The  two  cases  were  cho¬ 
sen  to  check  our  method  of  trajectory  calculation  for  both  high  and  low  collision  efficien¬ 
cies.  For  both  cases,  Langmuir’s  0  parameter  was  chosen  to  be  104.  This  is  given  by: 
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droplet  density 
air  density 

dynamic  viscosity  of  air 
cylinder  diameter 
freestream  speed. 


(9) 


0  is  a  nondimensional  impingement  parameter.  Large  values  of  0  imply  a  large  radius  of 
curvature  of  the  streamlines,  and  vice  versa.  Langmuir’s  A  parameter  was  36,0  in  the  first 
case  and  1.0  in  the  second.  A  is  given  by  the  expression: 


A  = 


9  MaC 


(10) 


where  rd  is  the  droplet  radius.  A  is  the  nondimensional  inertia  parameter.  It  is  the  ratio  of  the 
droplet’s  projectile  range  under  Stokes’  law  to  the  radius  of  the  cylinder.  If  A  is  small,  the 
droplets  tend  to  follow  the  streamlines,  and,  hence,  the  collision  efficiency  tends  to  be  low. 

In  these  cases,  as  in  all  the  experiments  considered  in  this  report,  the  droplets  were  intro¬ 
duced  into  the  airstream  five  chord  lengths  upstream  of  the  nose  of  the  target  with  a  Rey¬ 
nolds  number  Rcd  of  0.001 .  Ideally,  the  droplet  trajectory  integration  should  begin  infinitely 
far  upstream  with  the  droplets  having  the  same  velocity  as  the  air  (i.e.  Red  =  0),  but  for  com¬ 
putational  reasons  this  is  impractical.  Tests  indicate  that  the  trajectory  errors  caused  by  this 
imperfect  initial  condition  arc  smaller  than  the  numerical  integration  errors. 

The  parameters  chosen  for  the  two  cases  considered  arc  given  in  Table  L  Table  l  also  pre¬ 
sents  a  comparison  between  our  results  and  those  of  Langmuir  and  Blodgett  (1946).  The 
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Table  1.  Icing  on  a  cylinder,  present  calculations  (rows  2,  4,  5)  vs  Langmuir  and  Blodgett 
(rows  1,  3). 
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symbols  vx  and  vy  denote  the  droplet  impact  velocity  components  in  the  x-  and  ^-directions, 
respectively.  0m  denotes  the  maximum  angle  of  droplet  impingement  from  the  forward  stag¬ 
nation  point.  Our  model  results  for  vx,  vy,  £m,  /J0  and  0m  (given  in  rows  2  and  4)  compare  fa¬ 
vorably  with  those  of  Langmuir  and  Blodgett  (given  in  rows  1  and  3).  The  discrepancies, 
which  are  larger  for  the  case  with  the  smaller  inertia  parameter  A,  are  quite  acceptable,  if  one 
recognizes  that  the  Langmuir  and  Blodgett  data  should  not  be  viewed  as  an  absolute  stan¬ 
dard. 

Figure  I  shows  the  How  field  (indicated  by  velocity  vectors)  and  the  droplet  trajectories 
(indicated  by  solid  curves)  for  flow  about  the  cylinder  in  the  two  cases  just  considered.  It  is 
interesting  to  note  the  much  higher  curvature  of  the  trajectories  for  the  less  massive  drops 
and  the  larger  shadow  zone  between  the  grazing  trajectory  and  the  cylinder.  One  might  also 
speculate  on  the  collisions  that  could  occur  between  the  small  and  large  drops,  because  of  the 
way  the  smaller  ones  track  across  the  trajectories  of  the  large  ones.  It  is  hard  to  see,  however, 
how  such  collisions  might  have  any  significant  effect  on  the  icing. 

Figure  2  displays  the  collision  efficiency  for  the  two  cases  as  a  function  of  the  nondimen- 
sional  distances  along  the  cylinder  surface  from  the  stagnation  point  (that  is,  actual  distance 
L  divided  by  cylinder  diameter  C).  Negative  values  of  L/C  lie  below  the  stagnation  point, 
while  positive  values  lie  above  it.  Figure  2a  is  almost  a  cosine  curve,  a  result  that  would  occur 
if  the  droplet  trajectories  were  straight  lines.  The  slight  “kinks"  in  Figure  2b  are  artificial 
and  arise  from  the  numerical  spline  fitting  procedure.  The  overall  collision  efficiency  is  equal 
to  the  total  area  under  the  curves. 

Figure  3  is  similar  to  Figure  2,  except  that  the  abscissa  is  now  A VC,  where  X  is  the  projec¬ 
tion  of  the  arc  length  L  onto  the  v-axis.  The  effect  of  this  change  in  abscissa  is  to  “squeeze" 
the  curves  in  towards  the  origin.  This  squeezing  is  greatest  near  the  origin,  so  that  the  most 
apparent  effect  is  to  sharpen  the  peak  in  the  curves.  Although  it  is  somewhat  redundant  to 
present  collision  efficiencies  as  functions  of  X VC  and  L/C  for  a  cylinder,  the  difference  is 
more  meaningful  for  airfoils,  as  some  of  the  historical  papers  prefer  X VC  and  others  prefer 
L/C.  For  airfoils,  a  plot  of  collision  efficiency  vs  L/C  provides  more  resolution  near  the 
nose. 
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b.  C  -  15  cm ;  i  x  114.3  m  s  ;  rf/  7.0  pm . 

f  igure  3.  Collision  efficiency  vs  distance  along  chord. 


RESULTS  AND  DISCUSSION 

Comparing  results  with  and  without  the  history  term 

Rows  4  and  5  of  Table  1  present  results  of  numerical  experiments  run  respectively  without 
and  including  the  history  term.  All  the  other  conditions  of  the  experiment  are  identical.  The 
droplet  trajetories  calculated  with  the  history  term  are  less  influenced  by  the  rapid  changes  in 
the  airflow  just  before  the  collision,  and  so  they  tend  to  travel  along  straighter  paths  than 
those  whose  trajectories  ignore  the  history  term.  This  is  explained  by  the  fact  that  the  history 
term  acts  to  reduce  the  droplet  accelc  ation,  because  it  takes  into  account  the  finite  rate  at 
which  vorticity  is  shed  from  the  accelerating  droplet.  The  net  result  is  that  in  this  case,  ignor¬ 
ing  the  history  term  reduces  the  maximum  impingement  angle,  0m ,  by  3.3  °,  reduces  the  stag¬ 
nation  line  collection  efficiency,  /i0,  by  1.5%,  and  reduces  the  overall  collection  efficiency, 
Em,  by  1.7%  (about  10%  of  the  actual  value).  The  particular  case  used  to  study  the  influence 
of  the  history  term  was  chosen  to  give  a  large  effect.  In  most  cases,  the  effects  would  prob¬ 
ably  be  less  than  those  indicated  here,  suggesting  that  the  term  may  be  ignored  without 
severely  affecting  the  accuracy  of  the  calculations. 

Collision  efficiency  of  NACA  0015  airfoil  at  8°  attack  angle 

Thus  far,  the  computational  icing  experiments  have  been  limited  to  cylinders.  Let  us  now 
consider  the  case  of  a  NACA  0015  airfoil  at  an  attack  angle  of  8°.  The  chord  length  is  16.9 
cm,  the  freestream  speed  30.5  m  s'1,  and  the  droplet  diameter  20  ^m.  The  history  term  is  not 
included  in  the  computation.  Figure  4  illustrates  the  resulting  airflow  (indicated  by  velocity 
vectors)  and  droplet  trajectories  (indicated  by  solid  curves)  for  this  case.  It  should  be  kept  in 
mind  that  the  flow  region  depicted  in  the  figure  is  only  a  small  portion  of  the  total  flow  con¬ 
sidered.  In  addition,  the  coordinate  system  is  fixed  to  the  airfoil  so  that  the  flow  appears  to 
be  impinging  upwards  in  a  horizontally  oriented  airfoil.  In  fact,  the  entire  figure  should  be 
rotated  clockwise  by  8°. 

The  droplet  trajectories  clearly  indicate  the  asymmetry  of  the  impingement  above  and  be¬ 
low  the  stagnation  point,  when  the  airfoil  is  not  at  0°  attack.  This  is  generally  reflected  in  dif¬ 
ferent  icing  characteristics  above  and  below  the  stagnation  line.  Figure  5,  which  depicts  the 
local  collision  efficiency,  also  illustrates  this  asymmetry.  Negative  values  of  L/C  lie  below 
the  airfoil  nose  and  positive  values  above  it.  The  collision  efficiency  is  a  maximum  close  to 
(though  not  necessarily  at)  the  stagnation  line.  The  overall  collision  efficiency  for  this  case  is 
50,1%  and  the  maximum  is  74.4%  at  a  distance  L/C  =  -0.009  below  the  nose.  A  compari¬ 
son  of  Figure  5  has  been  made  with  the  results  of  Bragg  (private  communication),  who  has 
also  recently  investigated  this  problem  (see,  for  example,  Bragg  and  Gregorek  1981).  The  dif¬ 
ferences  between  the  two  sets  of  computed  results  are  generally  negligible,  in  the  sense  that 
experiments  would  not  likely  be  of  sufficient  accuracy  to  allow  one  to  choose  between  the 
two. 

Time-dependent  accretion  on  NACA  0015  airfoil  at  8°  attack  angle 

The  collision  efficiencies  plotted  in  Figure  5  are  those  for  the  airfoil  surface  itself  at  the 
onset  of  icing.  Once  a  significant  accretion  has  built  up  on  the  airfoil,  the  collision  efficien¬ 
cies  change,  and  this  change  affects  the  subsequent  development  of  the  accretion.  This  feed¬ 
back  process  between  the  accretion  and  the  airflow  and  droplet  trajectories  goes  on  continu¬ 
ously  in  nature.  We  decided  to  simulate  the  continuous  process  in  a  step-wise  fashion.  Thus, 
using  the  computed  initial  collision  efficiencies,  we  estimate  the  profile  of  the  ice  accretion 
after  a  finite,  but  small,  time  interval.  We  then  use  this  new  airfoil  profile  (including  the  al¬ 
ready  accreted  ice  layer)  to  determine  a  new  flow  field,  droplet  trajectories,  and  collision  ef¬ 
ficiencies.  After  that,  a  new  increment  to  the  accretion  is  once  again  calculated,  and  the  en¬ 
tire  process  is  repeated  for  as  long  a  total  period  as  desired.  Other  authors  (e.g.  Lozowski  et 
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al.  1979)  have  not  taken  this  feedback  process  into  account,  but  instead  have  used  the  initial 
collision  efficiencies  to  try  to  extrapolate  the  growth  of  the  accretion  over  substantial  periods 
of  time.  In  this  section  we  compare  these  two  procedures:  viz.  single-step  accretion  vs  multi- 
step  with  feedback.  We  also  make  a  comparison  with  an  experimental  accretion  grown  in  the 
NRC  high-speed  icing  tunnel. 

Figure  6  shows  the  airflow  and  the  droplet  trajectories  for  a  NACA  0015  airfoil  at  an  8° 
angle  of  attack,  but  the  conditions  are  somewhat  different  from  those  of  Figure  4.  In  the 
present  case,  the  chord  length  is  21 .3  cm,  thefreestream  speed  61ms"',  and  the  droplet  diam¬ 
eter  20  /un.  The  history  term  is  not  included  in  this  simulation.  The  solid  line  in  Figure  7 
shows  the  initial  local  collision  efficiency  0  vs  the  nondimensional  length  L/C  along  the  air¬ 
foil  surface.  Based  on  these  values  of  0,  and  assuming  a  cloud  liquid  water  content  of  0.4  g 
m-3  and  an  accretion  density  of  890  kg  m~\  the  accretion  growth  in  a  single  step  over  2.5 
minutes  was  calculated.  The  modified  airfoil  profile,  with  the  ice  accretion  attached,  was 
then  used  as  a  basis  for  calculating  a  new  airflow  and  new  droplet  trajectories.  From  these,  a 
new  determination  was  made  of  the  local  collision  efficiency  after  2.5  minutes  of  icing.  This 
is  indicated  as  the  dashed  curve  in  Figure  7  (L  now  being  measured  from  the  nose  of  the  ac¬ 
cretion  rather  than  from  the  nose  of  the  airfoil).  Although  the  differences  between  the  solid 
and  dashed  curves  are  not  striking,  it  is  clear  that  there  are  some.  In  particular,  £m  falls  with 
time  from  58.2%  to  56.5%  while  /30,  the  maximum  collision  efficiency,  actually  rises  from 
75.5%  to  78.9%.  Although  the  comparison  is  difficult  to  interpret  because  L/C  has  a  slightly 
different  meaning  in  each  case  (although  C  itself  remains  the  chord  length  of  the  basic  air¬ 
foil),  it  is  apparent  that  the  collision  efficiency  distribution  has  narrowed  and  become  more 
peaked  as  a  result  of  the  change  in  the  airflow  caused  by  the  first  2.5  minutes  of  accretion. 
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figure  6.  Trajectories  about  a  \.\C  \  00 15  an  foil  at  $  attack  am>le.  C  21.3  cm:  \  x  6/ 
m  s';  r({  --  10  fan . 
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/  inure  7.  c  ollision  efficiency  vs  length  alone  airfoil  surface  (SAC  \  0015  airfoil  at  <V  attack 
angle A  C  21.5  cm;  ~  61  m  s  \  (j  10  urn. 


The  decrease  in  the  collision  efficiency  that  occurs  with  time  on  the  lower  surface  near  the 
nose  is  illustrated  in  Figure  8,  which  shows  the  accretion  profiles  after  the  addition  of  two 
successive  2.5-minute  accretions.  This  two-step  accretion  shape  is  compared  in  Figure  9  with 
a  single-step  5-minute  accretion,  calculated  using  only  the  initial  collision  efficiencies.  The 
single-step  accretion  model  overestimates  the  growth  above  and  below  the  nost  and  underes¬ 
timates  it  at  the  nose  itself.  Although  the  differences  between  the  single  and  muhistep  models 
may  seem  relatively  minor  over  this  period  of  time,  they  will  be  much  more  significant  over 
longer  periods,  the  multistep  method  giving  much  more  realistic  results. 

The  shape  of  an  experimental  accretion  profile  made  under  similar  conditions  in  the  NRC 
high-speed  icing  tunnel  (Stallabrass  and  Lozowski  1978)  is  also  shown  in  Figure  9.  The  exper¬ 
imental  and  theoretical  results  are  not  perfectly  comparable  because  a  droplet  size  (approxi¬ 
mately  equal  to  the  medium  volume  diameter  of  the  tunnel  droplet  spectrum)  was  employed 
in  the  model.  The  general  agreement  is  quite  encouraging,  though  one  gets  the  impression 
that  the  model  accretion  occurs  too  low  on  the  airfoil  relative  to  the  experimental  one.  This 
discrepancy  rr  ay  be  the  result  of  a  bias  error  in  the  model.  On  the  other  hand,  it  may  have  to 
do  with  the  way  the  experimental  profiles  were  measured.  The  experimental  profiles  were  ob¬ 
tained  by  making  an  impression  in  plasticene  and  then  photographing  their  outline  against 
the  outline  of  the  airfoil.  Inaccurate  registration  of  the  airfoil  out'ine  and  that  of  the 
plasticene  mold  may  have  displaced  the  experimental  profile  upwards  from  where  it  should 
be.  Only  further  experiments,  with  an  improved  technique  for  measuring  the  experimental 
profile,  can  resolve  which  is  the  correct  explanation. 
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Figure  8.  Accretion  after  2.5  min  and  5  min  on  NACA  0015  airfoil  at  8°  attack  angle.  C  =  21.3  cm;  V ^ 
=  61  ms';  i(i  =  10  g.m;  LWC  =  0.4  g  m\ 
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Figure  9.  Accretion  after  5  min  in  a  single  step  ( solid )f  two  steps  (dashed),  and  experimentally  (hold) 
(NACA  00/5  airfoil  at  8°  attack  angle).  C  =  21.3  cm;  V  ^  =  61  ms';  x(j  =  10  ^m;  LWC  =  0.4  g  tn  \ 
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Time-dependent  accretion  on  NACA  0015  airfoil  at  0°  attack  angle 

In  this  section  we  demonstrate  that  the  multistep  accretion  process  can  be  continued  for  as 
many  as  five  steps,  after  each  of  which  a  new  potential  flow  and  new  droplet  trajectories  are 
calculated.  We  have  not  as  yet  attempted  to  increase  the  number  of  steps  beyond  five,  al¬ 
though  we  see  no  reason  why,  in  principle,  this  could  not  be  done.  Figure  10  illustrates  the 
airflow  and  droplet  trajectories  at  the  initial  time  before  ice  accretion  begins.  The  airfoil 
chord  length  is  21.3  cm,  the  freestream  velocity  91.5  m  s"1,  and  the  droplet  diameter  20  /on. 
To  be  rigorous,  the  history  term  was  included  in  the  trajectory  calculations  for  this  simula¬ 
tion,  although  generally  speaking  its  effect  may  not  be  large. 

Figure  1 1  shows  the  collision  efficiency  at  the  initial  time  and  after  3  minutes  of  ice  accre¬ 
tion.  The  local  collision  efficiency  increases  with  time  near  the  nose  and  diminishes  with  time 
farther  back  along  the  airfoil  surface.  Table  2  shows  that  the  overall  collision  efficiency  de¬ 
creases  with  time,  while  the  collision  efficiency  at  the  nose  increases  with  time. 

The  result  of  this  effect  on  the  accretion  itself  is  shown  in  Figure  12,  where  we  see  that  the 
accretion  tends  to  become  more  “pointed”  with  time,  and  that  the  growth  rate  at  the  nose  in 
the  model  increases  with  time.  This  result  seems  reasonable  inasmuch  as  the  effect  of  the  ac¬ 
cretion  is  to  decrease  the  local  radius  of  curvature  at  the  nose,  thereby  enhancing  the  colli¬ 
sion  efficiency  and  increasing  the  growth  rate.  Unfortunately,  no  time-dependent  experi¬ 
mental  growth  rate  measurements  are  available  to  confirm  this  result.  Figure  13  compares 
the  resulting  accretion  for  the  multistep  approach  with  that  obtained  using  a  single  5-minute 
step.  The  single-step  accretion  slightly  underestimates  the  growth  at  the  nose  and  greatly 
overestimates  it  farther  back. 

An  experimental  ice  accretion  profile  grown  under  similar  conditions  in  the  NRC  high¬ 
speed  icing  tunnel  is  included  in  Figure  13  for  comparison.  Although  the  agreement  is  good 
at  the  nose,  a  substantial  difference  occurs  farther  back.  We  suspect  that  this  is  due  to  the 
growth  of  low-density  feathery  rime  in  the  experimental  case.  Because  we  assume  an  ice 
density  of  890  kg  m"J  in  the  model,  feathery  rime  growth  is  not  taken  into  account. 


figure  10.  Trajectories  about  a  SAC  A  0015  airfoil  at  0  attack  angle.  (  21.  Jan;  1  ^ 

91.5  m  s  r(/  fO  /<////  L  WC  -  0.4  g  nt 
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COLLISION  EFFICIENCY  IN  PERCENT 


Figure  II.  Collision  efficiency  v.s  length  along  airfoil  surface  (NACA  0015  airfoil  at  0°  attack  angle). 
=  21.3  cm;  V ^  =  91.5  ms"1;  r(j  =  10  fim;  LWC  -  0.4  g  m~\ 


Table  2.  Collision  efficiencies 
as  a  function  of  time  for  the 
case  of  Figure  10. 


Time 

(min) 

(°’o) 

(°'o) 

0 

49.1 

81.8 

1 

47.7 

83.8 

46.1 

85.5 

3 

44.7 

86.3 

4 

43.3 

87.1 

17 


CONCLUSIONS  AND  RECOMMENDATIONS 


1.  The  principal  accomplishment  of  the  work  performed  under  the  present  contract  has 
been  the  development  and  testing  of  a  computational  simulation  model  of  rime  icing  on  arbi¬ 
trary  two-dimensional  airfoils.  The  computer  code  for  the  model  is  presented  in  Appendix  C. 
The  program  is  annotated  so  that  it  should  be  possible  for  scientists  elsewhere  to  run  the  pro¬ 
gram,  check  the  present  results,  and  develop  new  results  for  their  own  applications.  Should 
any  difficulties  be  encountered  in  the  implementation  of  the  model  program,  the  authors  will 
be  pleased  to  offer  their  advice  and  assistance. 

2.  Most  of  the  model  runs  described  in  the  present  report  have  been  performed  to  test  the 
accuracy  of  the  components  of  the  model.  The  potential  flow  field  was  tested  against  the 
known  analytic  solution  for  a  circular  cylinder,  and  was  found  to  behave  acceptably  (stream 
function  errors  of  0.1%  or  less,  which  give  rise  to  collision  efficiency  errors  of  less  than 
0.5%).  The  accuracy  of  the  droplet  trajectories  was  determined  by  comparing  the  model-pre¬ 
dicted  collision  efficiencies  with  those  computed  by  Langmuir  and  Blodgett  (1946).  Relative 
agreement  to  better  than  10%  was  found  even  in  the  worst  case.  Finally,  the  ice  accretion 
profiles  themselves  were  tested  against  two  experimental  accretions,  and,  while  the  agree¬ 
ment  was  not  exact,  it  was  encouraging  as  to  the  model’s  simulation  capabilities. 

3.  The  other  model  runs  presented  in  this  report  were  performed  either  to  test  the  impor¬ 
tance  of  the  history  term  in  the  droplet  equations  of  motion  or  to  compare  a  single-step  vs  a 
multistep  accretion  process.  Although  these  tests  were  not  exhaustive,  they  did  indicate  that 
omitting  the  history  term  did  not  have  a  dramatic  effect  on  the  results.  The  biggest  effect  of 
the  history  term  occurred  in  cases  with  low-  collision  efficiencies.  The  tests  also  showed  that 
the  accretion  profiles  predicted  by  the  single-step  and  multistep  processes  are  different,  and 
that  the  difference  increases  with  the  duration  of  the  accretion.  As  a  result  of  these  tests,  and 
because  in  principal  the  multistep  accretion  procedure  better  simulates  what  is  happening  in 
nature,  we  recommend  that  the  single-step  approach  be  used  only  for  brief  accretion  dura¬ 
tions.  Thus,  for  example,  a  single-step  model  might  be  quite  useful  for  helicopter  deicing  cal¬ 
culations.  On  the  other  hand,  the  multistep  method  would  be  preferable  for  simulations  of 
powerline  icing  where  the  duration  may  be  hours  or  days. 

4.  Within  the  scope  of  the  present  contract,  it  has  not  been  possible  to  use  the  model  to  in¬ 
vestigate  the  effects  of  various  parameters  on  the  shape  and  development  of  the  accretion. 
We  recommend,  however,  that  such  studies  be  undertaken  with  the  model.  Although  the 
model  is  presently  limited  to  simulating  rime  accretions,  there  are  many  questions  that  it  can 
be  used  to  investigate.  What  is  the  effect  of  airfoil  size  and  shape  on  the  accretion?  What 
would  be  the  effect  on  the  accretion  of  using  a  realistic  cloud  droplet  distribution  (see,  for  ex¬ 
ample,  Ackley  and  Templeton  1979)  rather  than  a  single  droplet  size?  What  parameters 
should  be  simulated  to  properly  model  ice  accretions  at  a  reduced  scale?  How  is  the  accretion 
changed  if  the  airfoil  attack  angle  and  the  airstream  speed  oscillate  as  they  would  for  a  heli¬ 
copter  rotor  blade?  Such  questions  and  many  others  could  and  should  be  profitably  con¬ 
sidered  using  the  present  icing  model. 
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APPENDIX  A:  SAMPLE  INPUT 

This  appendix  contains  a  sample  of  the  parameters  that  must  be  input  to  the 
program  along  with  examples  of  their  typical  values*  These  parameters  are  read 
by  the  program  from  input  device  4  (see,  for  example,  program  line  186). 
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APPENDIX  B:  SAMPLE  OUTPUT 

This  appendix  contains  a  sample  of  the  printed  output  produced  bs  the  pi  obtain 
starting  with  the  sample  input  values  given  in  Appendix  A.  This  output  is  wiiiien  to 
output  device  7.  The  trajectories  calculated  correspond  to  those  of  tiguie  I .  I  he  col¬ 
lision  efficiency  data  (on  the  last  page)  cot  respond  to  f  igures  3  and  5. 

FOR  EON.  SOLN.  IER  =  O 

THE  POTENTIAL  FLOW  LIFT  COEFFICIENT  IS  0.00000 
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APPENDIX  C:  PROGRAM  LISTING 


This  appendix  contains  the  program  listing  as  written  in  Fortran.  The  program  listing  has 
been  carefully  annotated.  However,  should  difficulties  be  encountered  in  attempting  to  run 
the  program  as  listed,  the  authors  are  prepared  to  offer  advice  and  assistance. 
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42 
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C 

C  WRITTEN  BY:  M.  OLESKIW  ON : 790526  LAST  MODI F I  ED : 80 1 228 
C 

C  CALCULATE  POTENTIAL  FLOW  ABOUT  AN  ARBITRARILY  SHAPED  AEROFOIL; 

C  CALCULATE  A  SERIES  OF  DROPLET  TRAJECTORIES  AND 

C  DETERMINE  THE  COLLISION  LOCATIONS;  FIND  THE  RESULTING  COLLISION 
C  EFFICIENCY  AND  ACCRETE  A  LAYER  OF  ICE. 

C  REPEAT  THE  PROCESS  FOR  A  PREDETERMINED  NUMBER  OF  STEPS. 

C 

C  INTERNAL  SUBROUTINES: 

C  COORDS:  CALCULATE  THE  UPPER  AND  LOWER  SFC .  COORDINATES 
C  OF  THE  AEROFOIL. 

C  POT  1  :  SOLVE  FOR  SFC.  VORTEX  DENSITY  ON  1  ELEMENT  AEROFOIL 
C  IN  POTENTIAL  FLOW,  GIVEN  COORDINATES  OF  AEROFOIL  SFC. 

C  STRMFN:  CALCULATE  STRE AMFUNCT ION  ON  A  GRID  ABOUT  AN  AEROFOIL 
C  SECTION  GIVEN  THE  SFC.  VORTICITY  DENSITY  ON  THE  AEROFOIL 

C  AND  PLOT  THE  FLOW  USING  VELOCITY  VECTORS. 

C  AIRPLT:  PLOTS  AEROFOIL  OUTLINE  WITHIN  WINDOW 

C  SFC .CALCULATE  Y  VALUES  AND  THE  LENGTH  FROM  THE  NOSE  ON  THE 

C  SFC.  OF  THE  AEROFOIL  BY  A  CUBIC  SPLINE  INTERPOLATION. 

C  SFCLEN: CALCULATES  THE  LENGTH  ALONG  A  SEGMENT  OF  A  CUBIC  SPLINE. 

C  CE -.CALCULATE  AND  PLOT  COLLISION  EFFICIENCY  OF  ARBITRARY 

C  AEROFOIL  BY  DETERMINING  A  SET  OF  IMPACTING  TRAJECTORIES. 

C  PLTSZ: DETERMINE  PARAMETERS  NECESSARY  FOR  SCALING  OF  A  PLOT 
C  AND  ITS  .AXES  . 

C  I C I NG : CALCULATE  AMOUNT  OF  ACCRETION  AND  DETERMINE  A  NEW  SET 
C  OF  AEROFOIL  SFC.  ELEMENT  ENDPOINTS  AFTER  DETERMINING  THE 

C  AEROFOIL  NOSE  LOCATION. 

C  GROWTHPLOTS  SUCCESSIVE  AEROFOIL  OUTLINES  WITHIN  VIEW  WINDOW. 

C  TRAJEC:  CALCULATES  TRAJECTORIES  OF  DROPLETS  IN  POTENTIAL  FLOW 
C  ABOUT  AN  AEROFOIL. 

C  ACCN : CALCULATES  RHS  OF  NON-DIMENSIONAL  EONS.  OF  MOTION. 

C  AIRVEL :CALCULATES  THE  AIR  VELOCITY  COMPONENTS  AT  A  GIVEN 
C  LOCATION. 

C  DRAG  CALCULATES  THE  REYNOLDS  NUMBER  AND  DRAG  COEFFICIENT 

C  OF  THE  DROPLET  AT  ANY  STEP  ALONG  ITS  TRAJECTORY. 

C  HIST iDETERMINES  VALUE  OF  INTEGRAL  IN  HISTORY  TERM. 

C  RKF4:  THE  RUNGE-KUTTA-FEHLBERG  4 TH  ORDER  ODE  INTEGRATION  TECHNIQUE 

C  RK4 :  THE  RUNGE-KUTTA  4TH  ORDER  ODE  INTEGRATION  TECHNIQUE. 

C  PC4 :  THE  PREDICTOR-CORRECTOR  4TH  ORDER  ODE  INTEGRATION  TECHNIQUE 
C  INTERNAL  FUNCTIONS: 

C  NSURF CALCULATES  THE  UNROTATED  X  VALUE  OF  A  POINT  ON  THE 
C  ACCRETED  AEROFOIL  SFC.  BASED  UPON  THE  COLLISION  EFFICIENCY, 

C  DIRECTION  OF  GROWTH,  AND  OLD  AEROFOIL  (ROTATED)  SFC.  POSITION 

C 

C  EXTERNAL  SUBROUTINES: 

C  I MSL :  (INTERNATIONAL  MATHEMATICAL  AND  SCIENTIFIC  LIBRARY) 

C  LEQT 1 F : SOLVES  SYSTEM  OF  EQNS . 

C  I CS I CU : CUB  I C  SPLINE  INTERPOLATION 

C  ZXGSN : GOLDEN  SECTION  SEARCH  METHOD  FOR  FINDING  FN.  MINIMUM. 

C 

C  SSPLIB.  (IBM  SUPPLIED  SCIENTIFIC  SUBROUTINE  LIBRARY) 

C  DELI  1  :  INCOMPLETE  ELLIPTIC  INTEGRAL  OF  THE  FIRST  KIND 

C  DELI2 : INCOMPLETE  ELLIPTIC  INTEGRAL  OF  THE  SECOND  KIND. 

C  DCE L  1  : COMPLETE  ELLIPTIC  INTEGRAL  OF  THE  FIRST  KIND. 

C  DCEL2 COMPLETE  ELLIPTIC  INTEGRAL  OF  THE  SECOND  KIND. 

C 

C  INPUT/OUTPUT  DEVICE  ASSIGNMENTS; 

C  3 : DAT  A  READ  BY  SUBPROGRAM  PLTSZ  TO  SCALE  PLOTS. 
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60  C  4  PROGRAM  INPUT  PARAMETERS  (DESCRIBED  BELOW). 

61  C  5  INPUT  CRT  OEVICF  FOR  CONTROL  OF  PROGRAM. 

G2  C  6  OUTPUT  CRT  DEVICE  FOR  MONITORING  OF  PROGRAM. 

63  C  7  OUTPUT  HARDCOPY  DEVICE  FOR  PRINTED  OUTPUT. 

61  C  9  OUTP'IT  FILE  FOR  STORAGE  OF  PLOT  DESCRIPTION  (CALCOMP  FORMAT). 

65  C 

66  C  PROGRAM  INPUT  PARAMETERS: 

67  C  TO  BE  READ  IN  FROM  INPUT  DEVICE  4.  EACH  GROUP  OF  PARAMETERS 

^8  C  IS  TO  BE  READ  FROM  THE  SAME  LINE  (CARD)  USING  THE  SPECIFIED 

69  C  FORMAT.  EACH  DATA  LINE  PRECEDED  BY  A  DESCRIPTIVE  REMINDER  LINE. 

70  C  SEE  EXAMPLE  FOR  DETAILS. 

7  1  C 

72  C  NE  F  =N0  OF  ELEMENT  FNDPTS .  ON  FRONT  HALF  OF  AEROFOIL  (14) 

73  C  NEB  =  N0 .  OF  ELEMENT  ENDPTS .  ON  BACK  HALF  OF  AEROFOIL 

74  C  (INCLUDES  THE  MIDPOINT  ENDPT .  (AT  THETA  =  90 ) )  (14) 

75  C  N I F  =  NO .  OF  SPLINE  ENDPTS/ELEMENT  ENDPT.  (FRONT  HALF)  (14) 

76  C 

77  C  ALPHA=ANGLE  OF  ATTACK  IN  DEGREES  (F6.0) 

78  C  TYPE =AEROFOI L  TYPE  (15) 

79  C  -1  ANALYTICAL  CYLlNOER 

80  C  O  NACA  RAZOR 

81  C  1  CYLINDER  (VORTEX  SHEETS) 

82  C  THICK-THICKNESS  OF  AEROFOIL  IN  PERCENT  (F6.0) 

83  C  X  M  I  N  = 

84  C  X  MA  X  -  VIEWPORT  SIZE  IN  X  (2F5.0) 

8  5  C  v  M I N  = 

86  C  V  MA  X  =  VIEWPORT  SIZE  IN  Y  (2F5.0) 

87  C  *7=  VELOCITY  VECTOR  GRID  SIZE  IN  X  (13) 

88  C  YZ=  VELOCITY  VECTOR  GRID  SIZE  IN  Y  (13) 

89  C  ANAL=0  ESTIMATE  SEGMENT  LENGTH  BY  NUMERICAL  APPROXIMATION.  (15) 

PC  C  1  : DETERMINE  SEGMENT  LENGTH  BY  ANALYTICAL  METHOD. 

q  i  r; 

92  C  PLTPAC=  PLOT  REDUCTION  OR  EXPANSION  FACTOR  FOR  ALL  PLOTS  (F7  2) 

93  C  TR.JPLA-PLOT  TRAJECTORIES  (0  OR  1)  (17) 

94  C  YOL  =  PLOT  THE  YO  VS  L  GRAPH  (0,  1.  OR  2 )  (2  PLOTS  AT  HALF  PAGE  SIZE)  (14) 

>5  C  CEL  =  PLQT  THE  CE  VS  L  GRAPH  (O,  1.  OR  2 )  (2  PLOTS  AT  HALF  PAGE  SIZE)  (14) 

96  C  CE  X "PLOT  THE  CE  VS  X  GRAPH  (O,  1.  OR  2 )  (2  PLOTS  AT  HALF  PAGE  SIZE)  (14) 

97  ''  ICEPLA-PLQT  AEROFOIL  AND  ICE  LAYERS  (0  OR  1)  (17) 

98  C  L  Y RMA X  - MA X .  NUMBER  OF  LAYERS  TO  ACCRETE  (17) 

99  C  CETOL =CRI TERION  (FOR  CHANGE  IN  CE  BETWEEN  ENDPTS.)  TO  DETERMINE 

I'M  C  WHETHER  OR  NOT  TO  CREATE  NEW  ENDPTS.  (F6.2) 

101  0  ICE  =FRACT I  ON  OF  CHORD  LENGTH  TO  BE  ACCRETED  PER  LAYER  ASSUMING 

102  C.  A  COLLISION  EFFICIENCY  OF  100%  (F7.2) 

103  C 

104  r  UINF-FREESTREAM  VELOCITY  (M/S)  (F6.0) 

105  C  C=CHORD  LENGTH  (M)  (F6.0) 

106  C  PINF  =  FOCESTREAM  PRESSURE  ( KPA  )  ( F6  O) 

107  C  TINF=FREESTREAM  TEMPERATURE  (C)  (F6.0) 

108  C  RD=DROPLET  RADIUS  (MICROMETERS)  (F6.0) 

109  C  A  1  = 

MO  C  B 1 =PARAMET ERS  FOR  PRED I CTOR -CORR ECTOR  FORMULAE  (2D10.0) 

111  C 

112  C  COSiDRAG  COEFFICIENT  FORMULATION:  (14) 

113  C  =0  ABRAHAM  ( 1970) 

M 4  C  = 1  RE  <  0  01:  STOKE'S  DRAG 

115  C  0.01  <  RE  <  5:  SARTOR  AND  ABBOTT  (1975) 

116  C  RE  >  5:  ABRAHAM  (1970) 

117  C  =2: LANGMUIR  AND  BLODGETT  (1945) 

118  C  TR JPRA  =  PR  I  NT  TRAJCTORY  INFO  (O  OR  1)  (17) 

119  C  PR  I  NT  I =N0 .  OF  STEPS  AT  WHICH  TO  PRINT  TRAJECTORY  INFO 

120  C  WITHIN  VIEWPORT.  (17) 

121  C  PLOT  I =N0  OF  STEPS  AT  WHICH  TO  PLOT  TRAJECTORY  WITHIN  VIEWPORT. 

122  C  (16) 

123  C  PR  I NTO=NO  OF  STEPS  AT  WHICH  TO  PRINT  TRAJECTORY  INFO 

124  C  OUTSIDE  VIEWPORT  (17) 

125  C  CRIT  ^CRITERION  (EXPRESSED  AS  %  OF  DROPLET  RADIUS)  USED 

126  C  TO  INDICATE  SUFFICIENTLY  CLOSE  DROPLET  APPROACH 

127  C  TO  DENOTE  COLLISION  ( F5  O) 

128  C  BETAO=ESTIMATED  LOCAL  COLLISION  EFFICIENCY  AT  STAGNATION  PT 

129  C  (F6.0) 

130  C 

131  C  NTR A JU  =  MANUAL  MODE:  NO  OF  TRAJECTORIES  PRINTED/PLOTTED  (17) 

132  C  -  AUTO  MODE :  NO  OF  TRAJECTORIES  DESIRED  ON  UPPER  SFC 

133  C  NTRAJL  =  AUTO  MODE  NO  OF  TRAJECTORIES  DESIRED  ON  LOWER  SFC  .  (17) 


134 

135 

136 

137 

138 

139 
1  40 
14  1 

142 

143 

144 

145 

146 
1  47 
1  48 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 
1G9 
170 
1  7  1 

172 

173 

174 

175 

176 

177 
1  78 
1  79 
180 
18  1 
1  82 

183 

184 

185 

186 

187 

188 

189 

190 
1  9  1 
192 
1  94 
194 
lot, 
1  °6 

107 

108 
199 

29Q 
201 
202 
209 
20  1 
o  5 

206 

207 


C  AT  =0 : ST  ART  TRAJECTORIES  AS  SPECIFIED  BY  INPUT  TERMINAL.  (13) 

C  1 : AUTOMATICALLY  DETERMINE  TRAJECTORY  STARTING  POINTS 

C  AFTER  FIRST  ONE  FOR  EACH  SFC . 

C  BOTH  =  0 :  SYMMETRICAL  AEROFOIL  AT  O  DEGREES  ATTACK  - 

C  CALCULATE  TRAJECTORIES  FOR  UPPER  SFC.  ONLY. 

C  1:  CALCULATE  TRAJECTORIES  FOR  BOTH  SFCS .  (15) 

C  EQN=0:  EON.  OF  MOTION  INCLUDES  TERMS  A  AND  B  (NO  INDUCED 
C  MASS  OR  BUOYANCY)  (14) 

C  1:  EON.  OF  MOTION  INCLUDES  TERMS  APRIME  AND  BPRIME 

C  2:  EON  OF  MOTION  INCLUDES  TERMS  APRIME.  BPRIME.  AND 

C  CPRIME  (HISTORY  TERM) 

C  PC= INTEGRATE  BY  RUNGE-KUTTA  (O)  OR  PRED I CTOR -CORRECTOR  (1) 

C  (AFTER  FIRST  3  INTERVALS)  OR  RUNGE -KUTT A - F EHLBE RG  (2)  (13) 

C  DTS -NON-D I M .  INITIAL  TIME  STEP  (F6.0) 

C  EPS=  FOR  ODE  INTEGRATION  TECHNIQUE  RUNGE -KUTTA-F EHLBERG : 

C  LOCAL  ERROR  DIVIDED  BY  LOCAL  STEP  SIZE.  ( D8 . 0 ) 

C  ACN -0 : DROPLET  INITIAL  VELOCITY  VECTOR  SLIGHTLY  GREATER  THAN 
C  THAT  OF  THE  AIR  AT  THAT  POINT.  (14) 

C  1 : DROPLET  INITIAL  ACCELERATION  WEIGHTED  BY  CHANGE  IN 

C  POTENTIAL  FLOW  FIELD. 

C 

C  XO=X  (UPSTREAM)  COORD  FOR  TRAJECTORY  STARTING  PTS.  (F10.0) 

C 

C  YO=Y  (OFF  AXIS)  COORDS  FOR  TRAJECTORY  STARTING  POINTS.  (F10.0) 
C  INPUT  ONE  FOR  EACH  SFC.  ( AUTO -TRAJECTORY  MODE),  OR  FOR  ALL 

C  THE  TRAJECTORIES  DESIRED  OTHERWISE. 

C 

5  FORMA T ( / , 3  1 4 ) 

10  FORMAT ( / , F6 . 0 , 15, F60.4F5. 0.213, 15) 

20  FORMAT ( / .F7.2.I7.3I4.2I7.F6.2.F7.2) 

30  FORMAT ( ' OTHE  ACCRETED  AREA  FOR  LAYER'. 13,'  IS'.FIO.S,/. 

THE  ACCUMULATED  ACCRETED  AREA  IS'.FIO.S) 

C 

DOUBLE  PRECISION  ALPHA.  XE (  101  )  ,  YE (  10 1 ) . LEN , YNNUR , XNNLR , 

PI , X . DF  LO A  T , LUC  lOI).LL(IOI) .XS.CETOL. ICE . YNNLR , ACCRU . ACCRL . 
XNP . YNP , XURTLP . XLRTLP . ACCR . ACCRT , INTU , INTL , XNNUR , 

. XU(  101  )  . YU( 101  I  , XL( 101  ) , YL( 101 ) . THICK , S30.C30. 

*LR( 101). YLR (101) .DSORT.XN. YN , BPARU ( 4 ) , BPARL ( 4 ) , CU ( 100.3). 

. CL(  100,3), XUR (  lOI  )  , YUR (  101 ) , ALPHAR , THETA, INTUP. INTLP 
C 

REAL  XMAX.XMIN, YMIN.YMAX.PLTFAC 

INTEGER  I  ,  J  .  TYPE  .  X  Z  .  YZ  .  r  JPL  A  .  NCOU  ,  NCOL  ,  IERU,  IERL.LYRM1  , 

. PLT . LAYER , LYRMAX , NCOL  1  L . YOL .ICEPLA, AT. BOTH, FAIL, ANAL . 

ATYPE , I ABS ,IU(51),IL(51). NEB. NEF.NIF.NIFP1.CEX.il. I J.NEU.NEL 

C 

COMMON  ALPHAR , P I / A E RO 1 /X E , YE/NOSE/XN, YN/FO I L/XUR , YUR , 

XLR . YLR/LG/LU, LL/LA /ANAL/ AER03/NC0U , NCOL /ROTP/C 30, S30 
. /GRID/XMIN, XMAX , YMIN. YMAX , XZ , YZ/SFCS/XU . YU , XL , YL 
,/SPL  I NE/CU  .  CL  /  AER04/NEU  ,  NEL/ENDS/ IU .  IL 
. /NNOSE/XNP . YNP , XURTLP , XLRTLP 
C 

C  INPUT  PARAMETERS: 

RE  AD ( 4 . 5 )NEF .NEB.NIF 

RE  AD ( 4 ,  10) ALPHA , TYPE , THICK . XM IN, XMAX . YMIN. YMAX. XZ.YZ, ANAL 
READ(4,20)PLTFAC.TRJPLA, YOL , CEL , CEX , ICEPLA , LYRMAX , CETOL , ICC 
P ! =3 .  14  159265358979324 
C  INITIALIZE  PARAMETERS 

ALPHAR  =  ALPHA  *P I / 1  802 

NCOU  =  NE  F  +  NEB 

NCOL  = NCOU 

*N-0  ODD 

VN - O  ODO 

ACCRT -O  no 

NIFP  1  NIF+ 1 

I  J  -  1 

ATYPE  “  I ABS ( TYPE  ) 

C 

C  CALCULATE  AEROFOIL  COORDS 

C  UPPER  AND  LOWER  COORDS  FOR  LEFT  HALF  OF  AEROFOIL 
DO  110  I  -  1  . NFF 
IU(  I  )  =  I J 
I  L  (  T  )  -  I  J 

DO  1 40  J~ 1  , NT  F  P 1 

THE  T  A  - P I / 2  DO  *  DF l OAT (  {  I  O  *NI F P 1 + J- 1  ) /DFLOAT (NEF*NIFP 1  ) 


1  T 


2oa 

CALL  COORDS (TYPE .THICK. THETA, X, YU{ Id) . YL( I J) ) 

209 

XU( I J)=X 

2  10 

XL ( I d ) SX 

2  1  1 

I J=IJ+1 

2  12 

140 

CONTINUE 

2  13 

1  10 

CONTINUE 

2  14 

C  UPPER  AND  LOWER  COORDS.  FOR  RIGHT  HALF  OF  AEROFOIL. 

2  15 

DO  150  1=1 .NEB 

2  16 

THETA=PI/2 .  DO*  (  1  . DO+DF  L0AT (1-1 ) /OF L0AT ( NEB- 1 ) ) 

217 

CALL  COORDS ( TYPE . THICK, THETA , X , YU( Id).YL(Id)) 

2  18 

XU( Id  )  =X 

2  19 

XL( Id)=X 

220 

I U (  NE F I  )  =  Id 

22  1 

I L ( NE  F+ I  )  =  Id 

222 

Id=Id+  1 

223 

150 

CONTINUE 

224 

NEU= I J- 1 

225 

NEL  =NEU 

226 

LAYERS  1 

227 

C 

228 

c  transform  these  coords  to  one  VECTOR  OF  LENGTH  NCOU+NCOL-1 

229 

c 

IN  CLOCKWISE  ORDER.  WITH  XE ( 1 ) =XE ( NCOL+NCOU- 1 )  -  THE  LEADING 

230 

ioo 

DO  102  I-1.NC0U 

231 

1 1  =  I U  (  I  ) 

232 

XE ( I  )  =  XU ( II  ) 

233 

YE (  I  )  =  YU  (  I  I  ) 

234 

102 

CONTINUE 

235 

NCOL 1 =NCOL -  1 

236 

00  104  1  =  1, NCOL  1 

237 

J  =  NC0U  +  NC0L -  I 

238 

1 1  =  I  L  (  I  ) 

239 

X E ( d  )  =  XL (  I  I  ) 

240 

vE(d )=YL(  II  ) 

24  1 

104 

CONTINUE 

242 

c 

243 

C  ROTATE  UPPER  &  LOWER  SFCS  BY  30  DEG  ABOUT  NOSE  IN  ORDER 

244 

c 

TO  FIT  CUBIC  SPLINES 

245 

c 

-  SEE  KENNEDY  &  MARSDEN  (1976) 

246 

c  DO 

NOT  ROTATE  IF  AEROFOIL  IS  A  CYLINDER. 

247 

I F ( ATYPE . EQ  1 )GOTO  200 

248 

S30  =  5 . D -  1 

249 

C3C=DSQRT ( 3 .DO) /2 .DO 

250 

GOTO  210 

25  1 

200 

S30=0 . DO 

252 

C30  = 1 • 00 

253 

2  10 

DO  320  1=1. NEU 

254 

XUR(  1  )  =  (  XU( I  )-xU(  1 ) )^C30+( YU( I )-YU(  1  )  )*S30 

255 

Y  UR ( I  )  =  ( YU(  I  )- YU(  1  ) ) *C30~ ( XU ( I ) -XU(  1 ) )  *S30 

256 

320 

CONTINUE 

257 

DO  330  I = 1 , NEL 

258 

<  L  R  (  I  )  =  (XL(  I  )-<L(  1  )  )  *C30- (YL(I)-YL(U ) *S30 

25Q 

YLRmMYLII)-YLm  )*C30M  XU  I  )  XL(  1  )  )  *  S30 

2G0 

3  30 

CONTINUE 

26  1 

C 

262 

6  SET  PARAMETERS  for  SPLINE  FITTING 

26  3 

I  F  <  AT  rPF  EQ  1  )  GO  TO  2  20 

264 

RPARUM)  =  i  DO 

5 

BPARUt  2  )  =1,  DO,  (  <UR(  2  )  -XUR(  1  )  )  *  (  (  YUR(  2  )  -  VUR(  1  )  )/(  XUR(  2  )  -  XUR 

2^6 

(  1  )  )  -  DSQRT  (  3  POM 

267 

BPARLH  3  )  =0  DC 

7  68 

R  R  A R U ( 4  ) -  0  DO 

2  r'  9 

BPARL  (  1  )  =  1  00 

BPARL( 2)6  DO.  (  <LR( 2  >  -XLR(  1  )  )*(  (YLR(2>-YLR(  1  )  )/<  XLR{ 2  )  - 

7  1 

*  L R 1  1  )  )  *PSQR! 1 3  DO) ) 

2  7  7 

RPARK3M6'  DO 

2  ?3 

B p A Rl  (41=0  DO 

2  7  ! 

GOTO 

7  7  5 

72° 

RPARUf  1 ) 2 0  00 

7  76 

RPARIM  2 ) =0  DO 

2  ’  7 

RPARUf  3  )-(.'  DO 

2  7  8 

B  P  A  R  U ( 4  ) - C  DO 

7  7  a 

P.PARl  (  1  )  =  0  DC 

7  86 

BPARL ( 2 ) - 0  DO 

78  i 

PPAR!  (  3  )  -  DO 

28  2 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 
29  7 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 
3  10 
3  1  1 
3  12 
3  13 
3  1  4 
3  15 
3  16 
3  1  7 
3  18 
3  19 
320 

32  1 

322 

323 

324 

325 

326 

327 

328 

329 
1 30 

33  1 

332 

333 
3  34 
335 
3  36 

33  7 
3  38 
3  39 

34  n 
3  4  1 
3  J  2 
9  4  3 
34  1 

3  15 
9  16 
<  1  7 
1  18 

4  19 
350 
16  1 
3C<  2 
7  6  } 
164 
.156 


BPARU4  )=0.D0 

C  FIT  CUBIC  SPLINES  TO  EACH  SFC . 

230  CALL  ICSICU(XUR, YUR , NEU , BPARU , CU , 100 , I ERU ) 

CALL  ICSICU1 XLR, YLR.NEL , BPARL ,CL ,  100. I  ERL) 

C 

C  CALCULATE  INTEGRAL  OF  UPPER  AND  LOWER  SFC.  PROFILES 
C  FIND  THE  LENGTHS  FROM  THE  NOSE  TO  VARIOUS  ENDPTS . 

LU (  1  )  =0  DO 
LL ( 1 ) =0 . DO 
INTU=0 . DO 
I NTL  =0 . DO 

DO  340  1=2, NEU 
XS -XUR ( I )-XUR( I-  1  ) 

CALL  SFCLEN( XS. LEN,CU( I  -  1 . 3 ) , CU( I  -  1 . 2 ) . CU ( I  -  1 .  1  )  ) 
LU( I  )  =  LU  (  1-1  )  +  LEN 

I  NTU=  I  NTU  +  (  (  (  CU(  1-1.3  )  *XS/4 .00*CU(  I-1,2)/3.DO)*XS 
+  CU( I  -  1 ,  1 )/2 . DO) *XS  +  YUR( I  -  1  )  )*XS 
340  CONTINUE 

DO  350  1=2. NEL 
XS  =  XLR( I ) -XLR(  I  -  1  ) 

CALL  SFCLEN(XS.LEN.CL(I-1.3).CL(I-1.2).CL(I-1.1)) 
LL <  I  )  =LL  <  I  -  1  MLEN 

I NT  L  = I  NT  L+ (  (  ( CL (  1-1 .  3  )  *XS/4 . DO+CL ( I  -  1 . 2 ) / 3  DO ) * XS 
+  CL(  1-1,1  )/2  00)*XS+YLR( 1-1 ) )*XS 
350  CONTINUE 

I F ( LAYER  EO  1 ) GOTO  400 
XNNUR  =  ( XN-XNP ) *C30+( YN-YNP ) *S30 
YNNUR  = ( YN-YNP ) *C30- ( XN-XNP )  *S30 
C  ACCRETION  AREA  FOR  UPPER  LAVER 

ACCRU  = I NTU - INTUP+  YNNUR*  XURTLP- XNNUR *  VNNUR/2  DO 
IF ( BOTH . EO  1 )GOfO  4 10 
ACCR  =  2  DO* ACCRU 
GOTO  420 

4 10  XNNLR- I  XN-XNP ) *C30- ( YN-YNP ) *$30 
YNNLRM YN-YNP ) *C30+( XN-XNP ) *S30 
C  ACCRETION  AREA  FOR  LOWFR  LAVER 

ACCRL= INTLP  I NT  L - YNNLR * XLRTLP ♦ XNNLR * YNNLR/2  DO 
ACCR  = ACCRU+ACCRl 
420  ACCRMACCRM  ACCR 
LYRM 1 = LAYER- 1 

WR  I  TE  (  6 . 30  It.  YRM  1  .  ACCR  .  ACCRT 
WRITE ( 7 . 30) LYRM 1 .ACCR. ACCRT 
400  I NTIJP  =  I  NTU 
INTLP- INTL 

C 

IF(  LAYER  GT.lVRMAx  AND  ICE  PL  A  EO  DC.OTO  121 
IF (LAYER  GT.LVRMAX  AND  I C F PL A . EO . 2 ) GOTO  130 
IF(  TYPE  GE  0  )  C  A  L  L  POM 
PLT  =  TRJPl  A  *  Y  0!  *  C  F  L  ♦  C  F  X+TCEPl A 
I F ( PL  T  EO  0 )GOTO  120 
IF ( LAYER  GT  1  )GOTO  125 
C 

c;  OPEN  PLOTTING 
CALL  PLOTS 
CALL  ME  T R I C (  1  ) 
f  ALL  ORGE  P ( 6  0.5-0.5.01 
CALL  r AC TOR< PL  T f  AC  ) 

C 

126  IF  (  TR.JPL  A  FO  •>  » <  ,0  I  0  12  1 

r  PLOT  AFPOfOIl.  OUTLINE  AND  VELOCITY  VICTORS 
1  3<>  FAIL  SI  RMf  Nf  T  vPE  ) 

12  1  >.  A  (  I  A  I  RP1  T  1  !  A  r  f  R  .  T  R,  IP  1  A  .  L  Y  RM  A  *  ) 

I  F  (  \  A  r  E  R  (  i  I  t  »  RM  A  *.  )  GO  T  D  3  7Q 

f  C  A!  CUI  A  I  F  DROPLET  I  R  A, if  0  TOR  I  fc  5 

1  ;■■■-!  ini.A-rp  F  0  UCAM  T  R  A  Jt  < *  (  T  V  P  F  ,  T  R  JPl  A  ,  1  H  I  f K  .  A  I  .  p  ( ■  U  m 

I  r  (  i  A  i  ER  if  1  ‘CAM  T  RAiJEK 
I  M  A  T  EO  O)r,0T(l  360 

C  CALCULATE  COM  I  SIGN  M  f  I  C  T  f.  NC  Y 

r  A  1  L  C  r  (  T  91  .Ml  /  F*  ,ri  TFAC.  1HILK.I  A.fkM 

r 

r  A'CRFTf  ICF  AND  F  [  NP  Nf  W  AEROFOIL  SHARI 
YAM  I  r  I  NG 1  ‘  F  T  01  .  I  r  F  ,  BO  T  H  .  F  A  I  t  I 


356 

357 

358 

359 

360 

36  1 
362 
353 

364 

365 

366 

367 

368 

369 

370 

37  1 
37? 

373 

374 

375 
3  76 
37  7 

37  3 
3  79 
380 
*8  1 

38  2 

38  3 

384 

385 
336 
987 
^88 

389 

390 
34  1 

39  ? 

'*94 
445 
1 9  r, 

14  7 
343 
i  ^  9 
!r0 

14)  1 

4'  ■  3 

■  1 

-1  6 
.  i '  - ; 


IFILAVER  E0  LYRMAX  AND. ICEPLA  EO.O)GOTO  360 
IF(FAIL.EO  1  )  GOTO  360 
LAYER  =  LAVER+  1 
GOTO  100 

c 

C  PLOT  SUCCESSIVE  AEROFOIL  SHAPES  ON  ONE  PLOT. 

370  CALL  GROWTH!  ICEPLA , LYRMAX , PLTFAC . TRJPLA  ) 

360  IF(PLT .NE .01CALL  PLOT < O O . . 999 ) 

STOP 

END 

C 

C 

SUBROUTINE  COORDS ( TYPE , T . THETA . X , YU. YL ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON . 790928  LAST  MOD  I F I  ED : 80 1 022 

C 

C  CALCULATE  THE  UPPER  AND  LOWER  SFC.  COORDINATES  OF  THE  AEROFOIL 

C 

DOUBLE  PRECISION  X . YU . YL . DSQRT . B . C . T , THETA , DCOS , E IM2 , E IM 1 , 

E I . DABS. A  tB . DSIN. XI , ETA . E . TC 

C 

INTEGER  TYPE . ATYPE . I ABS 

C 

COMMON  /JOUK 1/A , B . E I 

C 

C  IN  TYPE -AEROFOI L  TYPE 
C  IN  T-AEROFOIL  THICKNESS  T N  PERCENT 
C  IN  THE  T  A  =  ANGLE  FROM  X  AXIS 
C  OUT  X  =  X -COORD  OF  AEROFOIL  SFC. 

C  OUT  /U= 

C  OUT  YL  =  UPPER  &  LOWER  Y-COORDS  OF  AEROFOIL  SFC. 

C 

ATYPE  =  I ARS ( TYPE  ) 

IF  (ATYPE  EQ  UGOTO  101 

C 

C  CALCULATE  THE  SHAPE  OF  A  NAC A  AEROFOIL  MODIFIED  TO  HAVE  A  RAZOR-LIKE 
C  TRAILING  EDGE  BY  REMOVING  A  LINEARLY  INCREASING  AMOUNT 
C  FROM  X=0 . 3  TO  <-10 

C  RFF  GREGORY.  N.  8  P.G.  WILBY  (1973).  A  R.C.  PAPER  *1261 
C  ABBOTT.  I  H  8  A  E  VON  DOENHOFF  (1959),  THEORY  OF  WING  SECTIONS. 

5  TL  672  A 1 2  1959,  P113  8  32  1 

r 

C  CALCULATE  AEROFOIL  x  8  Y  COORDS  FOR  EACH  SFC 
X  =  (  1  DO- DCOS (  THETA )  )/2 .DO 

B  -O . 2969DO*DSORT ( X ) -0 .  126D0*X-0 . 35 16D0*  X  *X 
C -O . 2843D0*  «  * *30  1015D0*X*  *4 
i\J-T/0  2D2  *  (  E*C  l 

I  F  <  X  GT  O  31)0  )  f  U  -  YU  -  (  x  -  O  3D0  )  *  2  1 0  -  3  *  T  /O  7D0/0  2D2 
IF(v  GT  O  9399499999  )  YU "O  DO 
r  L  -  -  YU 
RF  TURN 


1  M 
1 '  .  4 
.11'. 

1  1  1 
1  1 
i  1  1 


( 

C  CALCULATE  THF  x  8  t  COORDS  OF  A  CYLINDER 
101  xU  1  DO  *  DCOS  (  IMF.  T  A  )  ) /2  DO 

YtJ ~DSQR  T  (  0  2500  (x  0  5D0  )  M  x  -O  .  5D0  )  ) 

I F ( *  GT  0 . 9999949999 ) YU^O  DO 

Y l.  -  YU 

RETURN 
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SOL  VF  FOR  MJRF  ACT  vO^Tfy  DENSITY  ON  1  FlEMFNT  AEROFOIL  IN  POTENTIAL 
FIOW.  GIVEN  COOPTJG  OF  AEROFOIL  SURFACE 

REF  KfNNFDY.  .J  I  K  0  J  MARSDEN  (1976).  CAN  AERO  8  SPACE  JOUR.. 
V?  2 .  *5.  P 2 4  3  2 66 

SUBROUTINE  LEOTIf  OF  MMSLDPITB  LINEAR  E  ON .  SOLN  .  FULL  STORAGE 
MODF .  5PACF  FCONOMI/ER  SOLN 


nOURI  f.  PRf  f  I  SION  <F  (  Id  )  .  YE  (  IOI  )  .  x  (  1  IO  1  )  .  VC  (  IOI  )  .  R(  101  )  . 
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465 
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470 
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477 

478 
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502 

503 


DATAN. DABS. DSIGN.DL0G,SI( 100) . C0( 100 ) . PI . CL . 

K(  101 . 101 ) , WK  A  R  E  A (  101 ) ,D(  100 ). XT . YT . DE , DELTA . 
DXC,DYC,B,A,R1S.R2S.R3S,T3.T1 ,T2, ALPHAR , DCOS . DSIN . DSQRT 

C 

INTEGER  N.N1.J.J1, I DGT , IER. I . NCOU . NCOU 1 .NCOL. Jd 

C 

COMMON  ALPHAR.PI/AER01/XE. YE/AER03/NCOU. NCOL/ AER02/XC.YC.R.D. SI , CO 

C 

10  F ORMAT (  ' OFOR  EON.  SOLN .  IER='.I3./. 

.  OTHE  POTENTIAL  FLOW  LIFT  COEFFICIENT  IS',F9.5) 

20  FORMAT* 'OCONTROL  PT  .  X  COORD.  Y  COORD.  SFC .  VEL.') 

30  FORMAT (  '  '  . I  6 . 5X . 2F 10 . 5 , F 1 1 . 5 ) 

C 

NCOU 1 =NCOU -  1 
N=NCOU 1+NCOL- 1 
N 1 =N+ 1 

C 

C  CALC.  ELEMENT  LENGTHS  (D)  AND  CONTROL  POINTS  (XC.YC) 

C  XE (  1 ) =XE ( 2 *NCO- 1 )  =  XE ( N 1 )  =  LE AD ING  PT .  X  COORD. 

DO  110  d=1,N 
d1=d+1 

XC(iJ)  =  (XE(J)+XE(J1  ))  *0.500 
YC(d)=(YE(d)+YE(d1 )) *0.500 

0 ( J  )  = DSQRT ( (XE(J1 ) -XE ( d ) )**2+(YE(J1 )-YE(J))**2) 

110  CONTINUE 

C 

C  FIND  TRAILING  POINT  COORDS.  XC ( N 1 ) , YC ( N 1 ) :  FIG. 5 
XT=XE ( NCOU )-(XC( NCOU 1 )+XC ( NCOU ) )*0. 500 
YT=YE ( NCOU )-( YC( NCOU 1 ) +  YC ( NCOU ) ) *0 . 500 
XC ( N 1  ) =XE ( NCOU  )  + 1  . D - 2  *  XT 
YC ( N 1  )  =  YE ( NCOU  )  + 1  .D-2*YT 
C 

C  FORM  MATRICES  K  AND  R:  EONS.  9  &  10 

C  DO  FOR  EACH  SFC  ELEMENT  J  (COLUMN  OF  K)  AND  ROW  OF  R 
DO  120  d= 1 ,N1 

R  (  d ) = YC ( d  ) * DCOS ( ALPHAR )-XC( J ) *DSIN( ALPHAR  ) 

IF( J.E0.N1  ) GO  TO  140 
d  1  = J*  1 
DE  =0 ( J ) 

C  CALCULATE  ANGLE  OF  ELEMENT  TO  X-AXIS 
C0( J  )  =  <  XE ( J  1  )  - XE ( d  )  )  /DE 
SI(d)  =  (YE(J1  )-YE(d))/0E 
0ELTA-DE/2  DO 
140  00  130  1=1, N1 

IF( J.EQ.N1 ) GO  TO  150 

C  FINO  DISTANCE  BETWEEN  CONTROL  PTS.  I  ANO  d 
DXC  =  XC (  I  )-XC(d) 

D YC  =  YC (  I  )  - YC ( d ) 

C  CALCULATE  COMPONENTS  OF  EON  9  AND  FIG  2 
B~DXC*COl d)+DYC*SI ( d  ) 

A-DYC*CO(d)-DXC*SI(d) 

R1S=A*A+(B*DELTA)*(B*DELTA  ) 

R2S=A*A+(B-DELTA)*(B -DELTA) 

R3S=A*A+B*BDELTA*DELTA 
IF ( R3S  LT  1  D -  30 ) GO  TO  160 
T  3  =  DA  T  AN( 2 . DO  *A*DELTA/R3S) 

GO  TO  170 

160  I F ( DABS ( A  )  LT  1  D-30)G0  TO  180 

T  3  -  DA  T  AN(  (  B  OF  L  T  A  )  /  A  )  -  D A  T  AN(  (  B  -  DE  L  T  A  )  /  A  ) 

GO  TO  170 

180  T  3  =  OS  I GN ( P I  , A  ) 

170  T  1  =  ( B*DELT A  )  *DLOG( R IS ) 

T2  =  (B'DELTA  >*DL0G(R2S) 

K(I,J)  =  (T1-T2*2  DO  *  A  *  T  3  -  4 . DO  *DE  L  T  A )/4 . 00/ PI 
GO  TO  130 

C  FOR  LAST  COLUMN  OF  K 
150  K  <  I . J  )  -  1 . DO 

130  CONTINUE 

120  CONTINUE 

I DG  T -  8 

CALL  L  EOT 1F(K,  1.N1.  101  .R. I DGT , WK AR  E  A .  IER) 

C  ON  OUTPUT.  THE  SOLN  IS  IN  R 

r 

C  CALCULATE  THE  LIFT  COEFFICIENT 
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CL  =0 . DO 

00  200  JJ  = 1 , N 

CL  =  CL -2 . DO*R( JJ ) *D( JJ ) 

200  CONTINUE 

WR I TE ( 6 , 10)  IER.CL 
WR ITE ( 7,10)  IER.CL 
WRITE(7, 20) 

C  OUTPUT  AEROFOIL  COORDS.  AND  $FC .  VELOCITY. 

00  210  JJ=1 ,N1 

WRITE(7,30)JJ, XC( JJ) . YC( JJ) ,R( JJ) 

210  CONTINUE 

RETURN 
END 
C 
C 

SUBROUTINE  STRMFN( TYPE ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:  800222  LAST  MODIFIED:  801229 

C 

C  CALCULATE  STRE AMFUNCT I ON  ON  A  GRID  ABOUT  AN  AEROFOIL  SECTION 
C  GIVEN  THE  SFC .  VORTICITY  DENSITY  ON  THE  AEROFOIL  AND  PLOT  THE 
C  FLOW  USING  VELOCITY  VECTORS. 

C  REF:  KENNEDY.  J.L.  &  D.F.  MARSDEN  (1976),  CAN.  AERO.  &  SPACE  JOUR. 

C  V  22,  *5 ,  PP  243-256 

C 

DOUBLE  PRECISION  ALPHAR, XE( 101 ) . YE ( 101 ) . XC( 101 ) . YC( 101 ) .GAMMA ( 101 ) 
,0(1 00 )  ,  S I  ( 1 00 ) . CO ( 100) , DBLE ,YUP1 ,YLP1 . YU, YL.ZZ.DEN.PJKE, PJKA , DD , 
PIO 
C 

REAL  PS  I ( 372  1  ) ,K(  101 ) , DELTA , P I , ALPHAS , SNGL . SCO ,SSI .X.Y.DXC.DYC, 

. XMIN.XMAX, YMIN, YMAX.B.A, R 1S.R2S.T3, A TAN, SIGN . T 1 . T2 , 

. R , ABS , LOG , FLOAT, SIN, CO S,R3S,DX,DY .DPX.DPY, XPAGE , YPAGE , 

XT IP, YTIP.XP1 , YP1 , YM1 , U . V , AHL . AHLEN , SORT 

C 

INTEGER  XZ,YZ,TYPE,J, I , M, XZ 1  ,YZ1 ,F , N.NCOU.NCOL , L ,  I  I 

C 

COMMON  ALPHAR. PI  D/AERO 1/XE , YE/ AER03/NC0U.NC0L/ AER02/XC, YC .GAMMA, D. 

. SI .CO 

. / GRID/ XMIN, XMAX .YMIN. YMAX.XZ, YZ/SRCH/DD, II 

C 

C  IN  TYPE  =  AEROFO I L  TYPE. 

C 

C  PLOT  BOUNDARIES 

CALL  NEWPEN( 1 ) 

CALL  OR  I G I N( 999 , 2 1 .0,  10.5,5.0,5.0) 

CALL  AX2EP(3.5,3,2,0,0.9) 

CALL  AXIS2(0. ,0. , 'X/C‘ .-3,21 . ,0. ,XMIN, ( XMAX -XMIN )/2 1 . ,3.5) 

CALL  AXIS2(21  .  ,0.  .  '  ' . -  1 . -  10 . 5 , 90 .  , 0 .  . 0 .  ,  1 . 75 ) 

CALL  AX2EP ( 1 .75.3.3,0, 1.1) 

CALL  AXI52(0.  ,0.  ,  'Y/C'  ,3.10.5.90.  , YMIN, ( YM AX -YMIN)/ 10. 5, -  1 .75) 

CALL  AXIS2(0.  .  10. 5,  '  ' ,  1  ,  -2  1  .  , 0 .  . XMIN, ( XMAX -XMIN )/2 1  .,3.5) 

C  CHANGE  TO  SECOND  PEN 
CALL  NEWPEN( 2 ) 

N=NCOU+NCOL - 2 
PI=SNGL( PID) 

C  ALPHAR=ANGLE  OF  ATTACK  IN  RADIANS 
ALPHAS=SNGL( ALPHAR ) 

C 

C  CALCULATE  STRMFN .  ON  GRID. 

DO  120  J=1  ,XZ 

X-XMIN  +  FLOAT ( J- 1  )/FLOAT( XZ- 1  ) * ( XMAX -XMIN) 

DO  130  1  =  1, YZ 

C  PSI  IS  STORED  IN  VECTOR  FORM  BY  COLUMNS. 

M=(  J-  1  )*YZM 

Y  =  YMAX -FLOAT ( I  -  1 )/FLOAT( YZ- 1  )  +  ( YMAX-YMIN) 

PSI ( M ) -0 . O 

IF ( TYPE  EQ  -1  )GOTO  135 
DO  140  L= 1 . N 

C  FIND  DISTANCE  BETWEEN  CONTROL  PT .  L  AND  GRID  PT .  I,J. 

DXC  =  X -  SNGL ( XC ( L ) ) 

DYC= Y -SNGL ( YC( L ) ) 

C  CALCULATE  COMPONENTS  OF  EON  9  AND  FIG  2 
DE L T A  =  SNGL ( D ( L  )  )/2 .0 
SCO=SNGL(CO( L ) ) 
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578 

SSI =SNGL( SI ( L ) ) 

579 

B=DXC*SCO+DYC*SSI 

580 

A=DYC*SCO~DXC*SSI 

58  1 

R1S=A*A+(B+DELTA)*(B+DELTA) 

582 

R2S=A*A+(B-DELTA)+(B-DELTA) 

583 

R3S=A*A+B*B-DELTA*DELTA 

584 

IF(R3S.LT. 1 .  E  - 30 ) GO  TO  160 

585 

T3=ATAN(2  0*A*DELTA/R3$) 

586 

GO  TO  170 

587 

160 

IF(ABS(A)  .LT.  1 . E -30 )G0  TO  180 

588 

T3=ATAN( (B+DELTA)/A)-ATAN( (B-DELTA 

589 

GO  TO  170 

590 

180 

T3=SIGN(PI . A ) 

591 

170 

T1=(B+DELTA)*L0G(R1S) 

592 

T2=(B-DELTA) *LOG( R2S ) 

593 

K ( L )  =  ( T 1-T2  +  2 .0*A*T3-4 . 0*DELTA )/4  J 

594 

PSI(M)=PSI(M) -SNGL ( GAMMA ( L ) ) *K(L ) 

595 

140 

CONTINUE 

596 

R=Y*COS( ALPHAS )-X*S IN( ALPHAS) 

597 

C  ASSURE  THAT  PSI  ON  AEROFOIL  *  0. 

598 

PSI(M)=PSI(M  )  +  R-SNGL  (  GAMMA  (  N-M  )  ) 

599 

GOTO  130 

600 

C 

601 

c  streamfn.  for  a  cylinder. 

602 

135 

DEN=  (  X  - 5  .  D -  1  )  *  *  2-*- Y  *  Y 

603 

IF (DEN. LT . 1 . D-70)GOTO  136 

604 

PSI (M)=Y-Y/4 . DO/ ( ( X -5 . D“ 1 )**2+Y*Y ) 

605 

136 

PSI ( M ) =0 . DO 

606 

130 

CONTINUE 

607 

120 

CONTINUE 

608 

C 

609 

XZ  1  =xz- 1 

6  10 

Y  Z  1  =  YZ -  1 

6  1  1 

F  “0 

612 

11  =  1 

6  13 

C 

6  14 

DO  200  U=2.XZ1 . 2 

6  15 

DX= ( XMAX -XMI N ) /FLOAT ( XZ 1 ) 

6  16 

X=XMIN+FLOAT( J- 1 ) +DX 

617 

DPX  =  2  1  . / F  LOA  T ( X  Z 1 ) 

618 

C  ARROWHEAD  TAIL  IN  FRAME  COORDS. 

6  19 

XPAGE=FLOAT( J  -  1 )*DPX 

620 

XR  1  =x+DX 

62  1 

C  CHECK  IF  CENTERED  DIFFERENCING  IS  OK 

622 

I F ( XP  t . LE . SNGL ( XE (  1  )  )  )GOTO  220 

623 

CALL  SFC(DBLE(XP1  ),YUP1 ,  1 ,0,ZZ) 

624 

CALL  SFC(DBLE(XP1  )  ,  YLP  1  ,0,0, ZZ) 

625 

F  =  F+  1 

626 

I F ( X . LE . SNGL ( XE (  1  ) )  )GOTO  220 

627 

CALL  SFC ( DBLE (X),YU, 1 ,0,ZZ) 

628 

CALL  SFC(OBLE(X), YL.O.O.ZZ) 

629 

F  =  F+  1 

630 

C 

63  1 

C  DO 

FOR  EACH  COLUMN  OF  ARROWHEAD  TAILS 

632 

220 

DO  210  1=2, YZ1 ,2 

633 

DPY=10.5/FL0AT( YZ1 ) 

634 

DY  =  ( YMAX- YMIN )/FLQAT(YZ1 ) 

635 

Y  =  YMAX -FLOAT ( 1*1 )*DY 

636 

C  ARROWHEAD  TAIL  IN  FRAME  COORDS. 

637 

YP AGE  = 1 0  3- F  LOAT ( 1-1  )*DPY 

638 

M=( J-  1  )*YZ+I 

639 

IF( F  .  LE  .  1  ) GOTO  230 

640 

YP  1  =Y-DV 

64  1 

YM 1 = Y  +DY 

642 

C  IS 

CENTERED  DIFFERENCING  OK? 

643 

IF ( YP  1  . GE . SNGL ( YU )  . OR . YM 1 . L E . SNGL ( YL 

644 

I F ( Y . GE . SNGL ( YU ) )GOTO  250 

645 

C  CHECK  FOR  LOCATION  WITHIN  AEROFOIL 

646 

I F ( Y  GT . SNGL ( YL ) )GOTO  210 

647 

C  FORWARD  DIFFERENCING  IN  Y 

64  8 

U=(PSI(M)-PSI(M+1 ) )/DY 

649 

GOTO  240 

650 

C  BACKWARD  DIFFERENCING  IN  Y 

65  1 

250 

U  = ( PS  I ( M- 1 )-PSI(M))/DY 
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690 
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696 
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700 
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703 

704 

705 

706 
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709 
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720 
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GOTO  240 

C  CENTERED  OIF  FERENC  I NG  IN  Y 
230  U"(PSI(M'1 )PSI(M+1 ) J/2.0/DY 

240  IFCF.EO  0>G0T0  260 

C  IS  CENTERED  DIFFERENCING  OK? 

I F ( Y . GE . SNGL ( YUP 1 ) . OR . Y . LE . SNGL ( YLP 1 ) )GOTO  260 
C  BACKWARD  DIFFERENCING  IN  X 

V=(PSI( ( J-2)*YZ+I  ) -PS I ( (0-1 )*YZ+I ))/DX 
GOTO  270 

C  CENTERED  DIFFERENCING  IN  X 

260  V  =  (  PS  I  (  (  J-2  )  *YZ+I  )-PSI  ( J*YZ+I  )  )/2  O/DX 

C  ARROWHEAD  TIP 
270  XTIP  =  XPAGE-HJ  +  DPX 

YTIP  =  YPAGE  +  V  +DPX 
AHL=SQRT(U*U+V*V ) 

C  ARROWHEAD  LENGTH 

AHLEN=0. 25*AHL*DPX 

CALL  AROHD ( XPAGE . YPAGE . XTIP. YTIP, AHLEN , O , 16) 

210  CONTINUE 

200  CONTINUE 

RETURN 
END 
C 
C 

SUBROUTINE  AIRPLT( LAYER . TRJPLA . LYRMAX ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON : 800607  LAST  MODIFIED:  801022 

C 

C  PLOTS  OUTLINE  OF  AEROFOIL  WITHIN  VIEW  WINDOW 

C 

DOUBLE  PRECISION  XU( 101 ) . YU( 101 ) . DO . XL( 101 ) , YL( 101 ) . 

. XE( 101 ) . YE( 101 ) 

C 

REAL  XMIN . XMAX . YMIN. YMAX , SNGL . XP . YP . XPT( 1 04 ) , 

. YPT( 104 ) , XPE ( 103 ) . YPE( 103 ) . XGR( 104 . 10) , YGR ( 104  , 10) . 

. XGRE (103,10). YGRE ( 103 . 10 ) . XPP . YPP 

C 

INTEGER  NCOU , NCOL , XZ , YZ,NCOB , IE , IP  , I  , J . NCOB 1,1,11, 

I T ( 10), LAYER, ITT , TRJPLA , I PB. LYRMAX . I TE ( 10 ) . I TTE . 

. IEL ,NEL,NEU,NELM2 

C 

COMMON  /GRID/ XMIN, XMAX .YMIN.YMAX.XZ . YZ/GROW/XGR . YGR , 

.XGRE .YGRE, ITE, IT/AER01/XE. YE/AERO 3/ NCOU, NCOL/ SR CH/DD, II 
. /SFCS/XU. YU, XL. YL/AER04/NEU.NEI 
C 

C  IN  LAYER= INDEX  OF  ACCRETION  LAYER 
C  IN  TR JPLA  =  PLOT  TRAJECTORIES  ( O  OR  1) 

C  IN  L  YRMAX= INDEX  OF  FINAL  ACCRETION  LAYER 

C 

NELM2 -NEL-2 
NCOB  =  NCOU*  NCOL -  1 
NCOB  1  =NCOB -  1 
IP^O 
I  E  =0 
C 

C  FOR  THE  UPPER  SFC . : 

DO  700  J= 1 , NEU 
X P  =SNGL ( XU( J  )  ) 

YP  *  SNGL ( YU( J )  ) 

IF(YP  GE . YMAXIGOTO  720 
t F ( XP  GE . XMAX )GOTO  730 
IP  =  IP+  1 
XPT(  I P  )  -  XP 
Y  P  T (  IP  )  =  Y  P 
700  CONTINUE 

GOTO  740 

720  I  F (  IP  GT  0 ) GOTO  750 
XPT  (  I  P *  1  )  =  XP 
YP  T (  IP*  1  )  =  YMA  X 

GOTO  760 

C  OUT  ALONG  THE  TOP  EDGE 

750  XPT<  IP* 1  W  XP  XPT(  IP  )  )/( YP  YPT(  IP  ))  M YMAX -YPT(IP))  +  XPT(IP) 
YPT (  IP*  1  )  =  YM A  X 
C  UPPER  RIGHT  CORNER 
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744 
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767 
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770 
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760  I P  = I  P  +  2 

XPT ( IP)  =  XMAX 
YPT< IP) =  YMAX 
GOTO  740 

C  OUT  ALONG  THE  RIGHT  EDGE 
730  XPT ( IP  +  1  )  =  XMAX 

YPT( I  P  +  1  )  =  (YP-YPT( IP) )/ ( XP-XPT ( I P ) ) * ( XMAX - XPT ( I P ) )+YPT(  IP) 

IP=IP+1 

C 

C  FOR  THE  LOWER  SFC  : 

740  I PB  =  I  P 

I  E  L  =0 

DO  800  J= 1 , NE LM2 
XP=SNGL (XL(NEL-J)  ) 

YP=SNGL( YL(NEL-J) ) 

I F ( XP . GE . XMAX  OR . YP . LE . YMIN )GOTO  820 
IF( J.EO.  1  )GOTO  830 

I F ( XPP . LE . XMAX . AND . YPP . GE . YMIN )GOTO  830 
IF(YPP.LE.YMIN) GOTO  840 
C  IN  ON  THE  RIGHT  EDGE 
IP=IP+ 1 
XPT (  IP  )  =  XMAX 

Y  PT  ( IP  )  =  ( YP-YPP )/( XP-XPP )  +  ( XMAX -XPP )  +  YPP 
GOTO  820 

C  IN  ON  THE  BOTTOM  EDGE 
840  XPT ( IP+1 ) = XMAX 

Y  PT (  IP  + 1  )  =  YM I N 

IP-IP+2 

XPT (  IP )  =  ( XP-XPP )/( YP - YPP )  M  YMIN-YPP )  +  XPP 
YPT ( I P ) = YMI N 
GOTO  820 

C  ADD  ANOTHER  POINT  WITHIN  WINDOW 
830  IP=IP+1 

X PT (  IP  )  =XP 
YPT  (  IP  )  = Y P 
820  XPP =  XP 

YPP  =  Y  P 

800  CONTINUE 

I F ( IP .NE . I PB ) GOTO  850 

IP=  IP  +  1 

XPT (  I P  )  *XMAX 

Y  P  T (  IP )=YMIN 
C 

C  ADD  PARAMETERS  NECESSARY  FOR  PLOTTING 
850  XPT ( IP+1  )  =XPT (  1  ) 

Y  PT (  IP+1  )  =  Y  PT (  1  ) 

X  P  T (  IP  +  2 ) =  XMI N 
YPT ( IP+2 )=YMIN 

DO  200  1=1 , NCOB 1 
X  P  =  SNGL ( X  E (  I  )  ) 

YP  =  SNGL ( YE (  I  ) ) 

IF(XP  GT.XMAX)GOTO  200 
IF(YP. GT.YMAX.OR.YP.LT. YMIN)GOTO  200 
I E  =  I E+  1 
XPE( IE  )  =  XP 
YP  E (  IE)  =  YP 
200  CONTINUE 

X  P  E (  IE+1  )  =  XM I N 
YPE(  IE+ 1  )  =  YM I N 
XPT (  I  P  +  3 )  =  ( XMAX -XMI N)/2 1  O 
XPE  (  IE  *2  )  MXMAX-XMIN)/2  1  -O 
YPT(  I  P  +  3  )  =  ( YMAX  YMIN)/ 10  5 
YPE(  IE  +  2  )-( YMAX-YMIN)/ 10  5 
I  T  (  LAYERHIPO 
ITT= IP+3 
C 

C  THESE  ARE  THE  AEROFOIL  OUTLINE  POINTS  TO  BE  PLOTTED  WITHIN  THE  WINDOW 
DO  400  I  =  1  .  I T  T 
XGR (  I  . LAYER  )  -XPT (  I  ) 

YGR (  I  ,LAYER)^YPT(  I  ) 

400  CONTINUE 

I T E ( LAYER  )  =  IE  +  2 
ITTE=IE+2 

C  THESE  ARE  THE  AEROFOIL  ELEMENT  FNDPTS  WITHIN  THE  WINDOW 


800 

801 

802 

803 

804 

805 

806 

807 

808 

809 

810 
81  1 
812 

813 

814 

815 

816 

817 

818 
8  19 
820 
82  1 
822 

823 

824 

825 

826 

827 

828 

829 

830 

83  1 

832 

833 

834 

835 

836 

837 

838 

839 

840 

84  1 

842 

843 

844 

845 

846 

847 

848 

849 

850 

85  1 

852 

853 

854 

855 

856 

857 

858 

859 

860 

86  1 
862 
863 
86  4 

865 

866 

867 

868 

869 

870 

871 

872 

873 


DO  450  1  =  1  .  ITTE 
XGRE( I ,LAYER)=XPE( I  ) 

YGRE( I ,LAYER)=YPE( I ) 

450  CONTINUE 

I F ( TR JPLA . EQ . 0 , OR . LAYER . GT . LYRMAX )GOTO  500 
CALL  NEWPEN ( 3 ) 

CALL  LINE(XPT,YPT, IP+1 , 1 ,0,0) 

CALL  L INEP ( 0 . 1 ) 

CALL  LINE(XPE.YPE, IE. 1 , -1 ,0) 

500  RETURN 
END 
C 
C 

DOUBLE  PRECISION  FUNCTION  NSURF(XROT) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON,  800905  LAST  MODIFIED:  801022 

C 

C  CALCULATES  THE  UNROTATED  X  VALUE  OF  A  POINT  ON  THE  ACCRETED  AEROFOIL 
C  SURFACE  BASED  UPON  JHE  COLLISION  EFFICIENCY,  DIRECTION  OF 
C  GROWTH,  AND  OLD  AEROFOIL  (ROTATED)  SFC .  POSITION. 

C 

DOUBLE  PRECISION  XUR( 101 ) . YUR< 101 ) . CU( 100 . 3 ) , XLR( 101 ) . YLR( 101 ) , 

. CL ( 10O , 3 ) , C30 , S30 . XROT , D . LENG . LEN , LU( 101 ) , LL ( 10 1 ) , 

. L( 5 1 ) . Y0(51 ) ,CEE( 50. 3) , XLRT , YLRT . XN , YN . DD . SLP , K . XLRN , YLRN , 

. DS I GN , DSORT . ICE , NSURF Y , CE 

C 

INTEGER  J , RUN , I , ICT, ICU, ICL , NEU , NE L , NE L 1 

C 

COMMON  /FOIL/XUR, YUR . XLR. Y LR/SPL I NE/CU , CL/R0TP/C3O , S30 
. /IND/NSURFY . ICE , I , J , RUN/LG/LU . L L/COL/L , YO, ICT. ICU, I  CL/E FF /CEE 
. /NOSE/XN. YN/AER04/NEU.NEL 
C 

C  IN  XROT  =  ROT ATED  X  POSITION  ON  LOWER  AEROFOIL  SFC. 

C 

10  F ORMAT (  ' OOUT  OF  BOUNDS  IN  SEARCHING  FOR  AEROFOIL', 

.'OR  CE  SPLINES  IN  NSURF') 

C 

I F ( J . LT .  1  )J=1 
RUN-RUN+ 1 
NEL 1 =  NE  L -  1 

C 

C  FIND  THE  APPROPRIATE  AEROFOIL  SPLINE  SEGMENT 
120  IF(XRQT,GT.XLR(J) )GOTQ  105 
J  =  J-  1 

IF(J.EOO)GOTO  600 
GOTO  120 

105  IF(XR0T.LE.XLR(J+1  ))GOTO  110 
J  =  J  +  1 

IF(J.LE.NEL1  ) GOT  0  105 
GOTO  600 

110  D=XROT-XLR( J ) 

C  FIND  LENGTH  ALONG  SFC.  FROM  NOSE  70  THIS  POINT. 

CALL  SFCLEN( D . LENG ,CL(J,3),CL(J,2),CL(J,  1  )  ) 

LEN=LL( J )+LENG 
C  ROTATED  COORDS. 

XLRT  =  XROT 

YLRT=YLR(U)+( (CL(J,3)*D+CL(U,2) )*D+CL( J .  1  )  )  *D 

C 

C  FIND  THE  APPROPRIATE  CE  VS  L  SPLINE  SEGMENT 
IF<  I  .  LT  .  1)1  =  1 

220  I F (  -LEN.GT . L(  I)  )GOTO  205 
1  =  1-1 

I F (  I  . EQ . O  )G0  T 0  200 
GOTO  220 

205  1 F< -LEN. LE . L( I + 1 ) )GOTO  210 

1  =  1+1 

IF<  I  LE  .  ICDGOTO  205 
GOTO  600 

C  CE  EQUALS  O  -  NEW  AND  OLD  SFCS  THE  SAME 
200  NSURF-XLRT ‘C30+YLRT  *  S30+  XN 

NSURF  v  =  -XLRT *  5  30  +  Y  L  R  T ♦ C30+  YN 
RETURN 

C 

C  CALCULATE  THE  COLLISION  EFFICIENCY 
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874 

875 

876 

877 

878 

879 

880 
881 
882 

883 

884 

885 

886 

887 

888 

889 

890 

891 

892 

893 

894 

895 

896 

897 

898 

899 

900 

901 

902 

903 

904 

905 

906 

907 

908 

909 

910 

91  1 

912 

913 

914 
9  15 

916 

917 

918 
9  19 
920 

92  1 

922 

923 

924 

925 

926 

927 

928 

929 

930 

93  1 

932 

933 

934 

935 

936 

937 

938 

939 

940 

94  1 

942 

943 

944 

945 

946 

947 


210  DD= -LEN-L ( I ) 

CE  =  ( 3 . D0*CEE ( I , 3)*DD+2  DO*CEE( 1,2)) *DD+CE E ( I . 1 ) 

C  FIND  AEROFOIL  SL^PE 

SLP= ( 3 . DO*CL ( U , 3  )  *D+2 . DO*CL( U,2))*D+CL(J. 1 ) 

K* -  1  DO/SLP 
C 

C  NEW  SURFACE  COORDS: 

XLRN=XLRT-DSIGN( DSQRT ( ICE* ICE*CE*CE/( 1  DO+K*K) ) ,K) 

YLRN=YLRT +K*  ( XLRN-XLRT ) 

NSURF  =  XLRN*C30+YLRN*S30+XN 
NSURFY  =  -XLRN*$30+ YLRN*C30+YN 
RETURN 

600  WRI TE ( 6 . 10 ) 

WR I T  E ( 7 ,  10) 

RETURN 

ENO 

C 

C 

SUBROUTINE  SFC ( X . Y , S , L . L EN ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:800623  LAST  MOD I F I  ED : 80 1 1 02 

C 

C  CALCULATES  Y  VALUES  AND  THE  LENGTH  FROM  THE  NOSE 
C  ON  THE  SFC.  OF  THE  AEROFOIL  BY  A  CUBIC  SPLINE  INTERPOLATION 

C 

DOUBLE  PRECISION  XN . YN . XUR ( 10 1 ) . YUR ( 10 1 ) . CU ( 100 . 3 ) , CL ( 100 . 3 ) , 
XLR( 101 ) , YLR( 101 ) ,XB, DELTA, DELTAP .DABS 
. . S30.C30.XR. YR.X, Y,LU( 101),LL( 10 1 ) . LEN . LENG . D 
C 

INTEGER  S , L , I , JU. JL.NEU1 , NEU . NEL 1 .NEL 

C 

COMMON  /NOSE/XN, YN/LG/LU . LL/FOIL/XUR . YUR . XLR . Y LR/SPL I NE /CU , CL 
. /R0TP/C3O, S30/AER04/NEU , NE  L/SRCH/D , I 
C 

C  IN  X  =  PO I  NT  AT  WHICH  Y  VALUE  IS  TO  BE  CALCULATED 
C  OUT  Y  =  S FC .  POSITION  ON  SPLINE 
C  IN  S  =  0 : LOWER  SFC. 

C  1  UPPER  SFC . 

C  IN  L  = 1 : F I NO  LENGTH  ALONG  AEROFOIL  SFC.  FROM  NOSE  TO  (X,Y) 

C  OUT  LEN  =  LENGTH  ALONG  AEROFOIL  SFC.  FROM  NOSE  TO  (X.  ) 

C 

10  FORMA T (  ' OOUT  OF  BOUNDS  ON  SEARCHING  FOR  SFC.  POSITION  '  , 

' IN  ROUTINE  SFC '  ) 

C 

JU  =  1 
JL  =  1 

C  ROTATED  X  COORD. 

XR  =  ( X - XN  )  *  C30 
IF(S . £00 )GOTO  150 

C 

C  FOR  THE  UPPER  SFC. 

NEU1=NEU-t 

I F ( XR . GT . 0 . DO )GOTO  120 
I F ( XR . LT  . O  DO  )GOTO  600 
Y  =  YN 

LEN=0  DO 
RETURN 

C  FIND  THE  APPROPRIATE  SPLINE  SEGMENT. 

120  IF(XR,GT  .  XUR  ( I ) )GOTO  105 

1  =  1-1 

I F ( I  EO . 0  )GOTO  600 
GOTO  120 

105  IF ( XR  LE . XUR( 1+1 ) )GOTO  110 
1  =  1+1 

I F ( I . LE . NEU 1 )GOTO  105 
GOTO  600 

1  10  D  =  XR -XUR (  I  ) 

C  ROTATED  Y  COORD 

YR  =  ( (CU(I  , 3  )  *D  +  CU (  1.2)  )*D  +  CU( 1,1)  )*D  +  YUR( I  ) 

XB  =  XR  *  C30- YR  *  S30+XN 
DELT  A  =  X-XB 

I  F( DABS (DELTA  )  . LE  ID- 10)G0T0  400 

DE L T AP  = -C30+S30* (  ( 3  D0*CU( I  , 3 ) *D  +  2 . DO*CU<  I , 2 ) )'D  +  CU( 1.1)) 
XR-XR-DELTA/DELTAP 


948 

UU  =  JU  +  1 

949 

GOTO  120 

950 

C 

95  1 

C 

UNROTATED  Y  COORD. 

952 

400 

Y  = YR*C30+ YN+XR  *  S30 

953 

IF(L.EO.O)GOTO  300 

954 

C 

FIND  THE  SEGMENT  LENGTH. 

955 

CALL  SFCLENt  D , LENG,CU( I , 3 ) . CU< I . 2 ) . CU(  I . 1 ) ) 

956 

LEN=LU(  I  )  +  LENG 

957 

GOTO  300 

958 

C 

959 

C 

FOR 

THE  LOWER  SFC . : 

960 

150 

NE  L 1 =NE  L -  1 

96  1 

c 

FIND  THE  APPROPRIATE  SFC .  SPLINE  SEGMENT. 

962 

220 

IF(XR.GT.XLR( I ) )GOTO  205 

963 

1  =  1-1 

964 

I F ( I . EQ . 0 )GOTO  600 

965 

GOTO  220 

966 

205 

IF(XR.LE  .XLR(  1+1  ) )GOTO  210 

967 

1  =  1  +  1 

968 

I F ( I .  LE . NEL 1 ) GOTO  205 

969 

GOTO  600 

970 

2  10 

D  =  XR-XLR(  I  ) 

97  1 

C 

ROTATED  Y  COORD. 

972 

YR=( (CL( I , 3)  *0  +  CL ( 1.2) )+D+CL( 1.1)) *D+ Y  LR( I ) 

973 

XB-XR*C30+YR  *S30+XN 

974 

DE  L  T  A  =  X-XB 

975 

IF(DABS( DELTA)  LE. 1 . D-10)GOTO  500 

976 

DEL  TAP = -C30-S30* ( ( 3 . DO*CL( I  , 3  )  *D  +  2  . DO*CL( 1.2) )*D+CL( 1,1)) 

977 

XR=XR-DELTA/DELTAP 

978 

JL  =  JL+ 1 

979 

GOTO  220 

980 

C 

98  1 

C 

UNROTATED  Y  COORD 

982 

500 

Y  =  -  X  R  +  S  30+  YR  *C30+  Y N 

983 

I F ( L  EO . 0 )GOTO  300 

984 

c 

FIND  THE  SEGMENT  LENGTH. 

985 

CALL  SFCLEN( D, LENG. CL( I,3),CL(I.2),CL(I,1)) 

986 

LEN  =  LL ( I )  +  LENG 

987 

300 

RETURN 

988 

C 

989 

600 

WR  I  T E ( 6 ,  10) 

990 

WRITE! 7 . 10) 

99  1 

RETURN 

992 

END 

993 

c 

094 

c 

995 

SUBROUTINE  SF CL E N ( D . L , A , B , C ) 

996 

c 

997 

c 

WRITTEN  BY:  M.  OLESKIW  0N:8OO525  LAST  MOD  I F I  ED : 800902 

998 

c 

999 

c 

CALCULATES  THE  LENGTH  ALONG  A  SEGMENT  OF  THE  CUBIC  SPLINE  FIT 

1000 

c 

AEROFOIL  SFC. 

10C1 

c 

1002 

c 

REF  : 

DOUG  S .  PHILLIPS  ( 1980) 

1003 

c 

1004 

DOUBLE  PRECISION  I I , NU , E . F . DSQR1 . DELT A , G , A , B , C , D , L , 

1005 

T1.T2.T3.T4. NU 1 . ANU 1 . DABS , NUO . ANUO , K , E2 , F2 . E3 . F3 , E02 . F02 . 

1006 

E03. F03. XO. X 1 ,CK . FO. EO. F 1 , E 1 . YP . D I STP . D I  ST . DF LOAT , Y . 

1007 

DLOG.DSIGN 

1008 

c 

1009 

INTEGER  IER, I .ANAL 

1010 

c 

101  1 

COMMON  /LA/ANAL 

1012 

c 

1013 

c 

IN 

D  =  ROT ATED  X  COORDINATE  OF  POINT  FROM  BEGINNING  OF  SEGMENT 

1014 

c 

OF  INTEREST  TO  WHICH  THE  LENGTH  IS  TO  BE  FOUND. 

1015 

c 

OUT 

L  =  Sf GMENT  LENGTH 

1016 

c 

IN 

A  = 

1017 

c 

IN 

B  " 

1018 

c 

IN 

C  =  SPLINE  PARAMETERS  FOR  SECTION  OF  INTEREST 

1019 

c 

44 


1020 
102  1 
1022 

1023 

1024 

1025 

1026 

1027 

1028 

1029 

1030 

1031 

1032 

1033 

1034 

1035 

1036 

1037 

1038 

1039 

1040 
104  1 

1042 

1043 

1044 

1045 

1046 

1047 

1048 

1049 

1050 

1051 

1052 

1053 

1054 

1055 

1056 

1057 

1058 

1059 

1060 
106  1 
1062 

1063 

1064 

1065 

1066 

1067 

1068 

1069 

1070 

107  1 

1072 

1073 

1074 

1075 

1076 

1077 
10"/ 8 

1079 

1080 

108  1 
1082 

1083 

1084 

1085 

1086 

1087 

1088 

1089 

1090 
109  1 

1092 

1093 


I I ( NU , E , F  )  =  NU/3  .  D0*DSQRT  (  1  .DO+(DELTA+NU*NU)**2  )* 

. ( 1 . D0+2 . D0*DELTA*G*G/ ( 1 . D0+NU*NU*G*G ) ) 

. +( ( 1 . D0+DELTA*G*G ) *F -2 . D0*DELTA*G+G*E )/3.D0/G**3 
C 

I F ( ANAL . EQ • 0 ) GOTO  200 
IF( A .NE . O . DO ) GOTO  100 
I F ( B . NE . O . 00 ) GOTO  110 

C 

C  A  AND  B  EQUAL  TO  O 

L=D*DSQRT( 1 .DO+C+C) 

RETURN 

C 

C  A  EQUAL  0.  B  NOT  EQUAL  O 

1 10  T  1  =  (  2  . DO*B  *  D  +  C  ) *DSQRT ( 1 . D0+ ( 2 . DO*B*D+C ) * *2 ) 

T2=C*DSQRT ( 1 .DO+C*C) 

T  3  =  DLOG ( ( 2 . DO*B  *D+C )+DSQRT ( 1  . D0+ ( 2 . DO*B*D+C ) * *2 ) ) 

T  4  =DLOG( C+DSQRT (  1 .DO+C *C) ) 

L=(T 1-T2+T3~T4)/4 .  DO/B 
RETURN 

C 

C  A  NOT  EQUAL  0 

100  NU 1 =D5QRT ( 3 . DO*DABS ( A ) ) * ( D+B/3 . DO/ A ) 

ANU 1 =DABS ( NU 1 ) 

NU0  =  B/3  D0/A+0SQRTO  .DO*DABS(  A)  ) 

ANUO=DABS ( NUO ) 

DELTA-(C'-B'fB/3.  DO/ A  )  *DS IGN(  1  DO,  A) 

G  = 1 . DO/ ( 1 . DO+DELTA*DELTA )**0. 25DO 
K=D$QRT ( 5 . D- 1 -DELTA  +G  +  G/ 2 . DO ) 

E2=0 . DO 
F2=0 . DO 
E02  =0 . DO 
F02  =  0 . DO 

X0=2  DO*G*ANUO/( 1 . D0'ANU0+ANU0*G*G ) 

X 1 =2 . DO*G* ANU 1 / (  1  . DO- ANU 1  * ANU 1 *G*G ) 

CK  =  DSQRT (  1  . DO-K  *K ) 

I F ( ANU 1 . EQ . 1 . DO/ G ) GOTO  120 
I F ( ANU 1  GT . 1 .DO/G)GOTO  130 
C 

C  ZETA  LESS  THAN  PI/2 

CALL  DELI  1 ( F 1 , X 1 , CK ) 

CALL  DELI2(E1 ,X1  ,CK,  1  . DO , CK *CK ) 

GOTO  140 
C 

C  ZETA  GREATER  THAN  PI/2 
130  CALL  DELI  1 ( F2 , -X 1 , CK ) 

CALL  DELI2(E2,-X1,CK, 1  DO , CK*CK ) 

C 

C  ZETA  EQUALS  PI/2 

120  CALL  DCEL 1 ( F 3 , K ,  I ER ) 

CALL  DCEL2(E3,K, 1 . DO , CK*CK , I ER ) 

F  1  =  2 . DO*  F3- F  2 
E  1  =2 . DO*E3-E2 

140  IF( ANUO.EQ. 1  D0/G)G0T0  150 
IF( ANUO.GT . 1 . DO/ G )G0T0  160 
C 

C  ZETA  LESS  THAN  PI/2 

CALL  DELI i(FO.XO.CK) 

CALL  DELI2(E0,X0,CK,  1  . DO , CK*CK  ) 

GOTO  170 
C 

C  ZETA  GREATER  THAN  PI/2 
160  CALL  DELI  1 ( F02, -XO,CK) 

CALL  DEL  I  2  I E02 , -XO.CK.  1  DO.CK^CK ) 

C 

C  ZETA  EQUALS  PI/2 

150  CALL  DCELHF03.K,  IER) 

CALL  DCE L2 ( E03 , K ,  1  DO ,CK*CK, IER) 

F0  =  2 . DO*  F03-F02 
E0=  2 . DO*  E03 -E02 

170  L  =  (DSIGN(  1  . DO. NU 1  ) *  I  I ( ANU 1 . E 1 . F 1 )  DSIGN(  1  DO. NUO)* I  I ( ANUO. EO. FO) ) 
. /DSQRT ( 3 . DO*DABS( A ) > 

RETURN 

C 

C  NON-ANALYTICAL  (APPROXIMATE)  SFC 


LENGTH  DETERMINATION 


1094 

1095 

1096 

1097 

1098 

1099 
1  100 
1  101 
1  102 
1  103 
1  104 
1  105 
1  106 
1  107 
1  108 
1  109 
1  1  10 
1111 
1112 

1113 

1114 

1115 
1  1  16 
1117 
1  1  18 
1  1  19 
1  120 
112  1 
1  122 
1  123 
1  124 
1  125 
1  126 
1  127 
1  128 
1  129 
1  130 
113  1 
1  132 
1  133 
1  134 
1  135 
1  136 
1  1  37 
1  138 
1  139 


200  L  =0 . DO 
YP  =  0 . DO 
D I 5TP=0 . DO 

DO  210  1=1.25 
DIST=D*DFLOAT( I )/25.DO 
Y=( ( A*DIST+B)*DIST+C)*DIST 
L=L+DSQRT ( (DIST-DISTP)**2+(Y-YP)**2) 

YP  =  Y 

D I STP  =  D I  ST 
210  CONTINUE 

RETURN 
END 

C 

C 

SUBROUTINE  CE ( YOL . CEL , CEX . PLTFAC . THICK , LAYER) 

C 

C  WRITTEN  BY:  M .  OLESKIW  0N:800622  LAST  MODI F I  ED : 80 1 227 
C 

C  CALCULATE  ANO  PLOT  COLLISION  EFFICIENCY  OF  ARBITRARY  AEROFOIL 
C  GIVEN  A  SET  OF  IMPACTING  TRAJECTORIES 

C 

DOUBLE  PRECISION  D , L ( 5 1 ) , Y0( 5 1 ) , BPAR ( 4 ) . CEE ( 50 . 3 ) . THICK , 

PN .P.DIST.SLP. SSLP, DABS , CET . ALPHAR , DCOS , CEMAX , PNI , ZZ.DCOS, 

LU  (  101  ),LL(  101 ) .XU(  101 ),XL(  101 ),YU(  101 ),YL(  101 ), Y . DBLE 

C 

REAL  LPMIN. YOPMIN . LRG , YORG . SNGL , FACT ( 2) , LP( 103) . FLOAT. 

. YOP( 1 03 ) .CEP1203) ,XPAR(4, 10).YPAR(4, 10 ) , LS ( 53 ) , YOS ( 53 ) , 

.PLTFAC , XP( 203) . XPM IN , CEPMIN , X . XLF .XRG.COS 

C 

INTEGER  CEL . F , I , ICT.IER.IRX.IRY.PX.PY. YOL . I CU . I CL . J , CEX  f II , 

. KK . KL , KU . L A Y  E  R , NEU , NE  L , CO , I J 

C 

COMMON  ALPHAR/COL/L. YO, ICT. ICU, I CL/E F F/CE E/PLTPRM/XPAR , YPAR 
. /CEM/P/LG/LU. LL/SFCS/XU, YU, XL . YL/SRCH/D , I J/AER04/NEU , NEL 
C 

C  IN  YOL  =PLOT  YO  VS  L  GRAPH( 0  OR  1) 

C  IN  CE  L  =  PLOT  CE  VS  L  GRAPH  (O  OR  1). 

C  IN  CEX  =  PLOT  CE  VS  X  GRAPH  (0  OR  1). 

C  IN  PLTF  AC  =  F ACTOR  FOR  SCALING  ALL  PLOTS. 

C  IN  THICK=AER0F0IL  THICKNESS  IN  %. 

C  IN  LAYERS  I NO EX  OF  ACCRETION  LAYER. 

C 

10  FORMAT (' -BETAO  (MAX  LOCAL  CE  )  IS',F7.1.'%  AT  A  DISTANCE  OF '  , 

. F 1 0 • 3 , '  FROM  THE  NOSE',/.'  THE  TOTAL  COLLISION  EFFICIENCY  IS', 
. F7 . 1 . ) 

20  FORMAT( ' -FAILURE  TO  CONVERGE  UPON  MAX  CE ' ) 


1  140 

114  1 
1  142 
1  143 
1  144 
1  145 
1  1  46 
1  1  M 
1  148 
1149 
1  150 

115  1 
1  152 
1  153 
1  154 
1  155 
1  156 
1  157 
1  158 
1  1  59 
1  160 
M6  1 
1  162 
1  163 
1  164 
1  165 
1  166 
1167 


C 

F  ACT (  1  )  =  1 .0 
F  ACT ( 2 ) =0 . 7 

C  CUBIC  SPLINE  END  PARAMETERS 
BPAR (  1  )  =  1  . DO 

BPAR ( 2 ) =6 . DO* <  Y0( 2 )-Y0<  1 )  )/( L( 2  )-L<  1  )  )*  *2 
BPAR ( 3  )  =  1  DO 

BPAR (4 )  =  -6.D0M Y0(  ICT)-Y0( ICT- 1  )  )/(L<  ICT)-L(  ICT  -  1  )  )*+2 
C  CREATE  SINGLE  PRECISION  VERSIONS  OF  L  AND  YO  IN  LS  AND  YOS 
DO  130  1  =  1  , ICT 
L S ( I  )  =  SNGL ( L ( I  ) ) 

130  CONTINUE 

DO  140  I  =  1  . ICT 
YOS ( I ) = SNGL ( Y0( I ) ) 

140  CONTINUE 

C  FIT  CUBIC  SPLINE  TO  YO  VS  L  CURVE 

CALL  ICS ICU ( L . YO. ICT .BPAR .CEE .50. IER ) 

C 

C  FIND  BETAO  (MAX  VALUE  OF  LOCAL  CE) 

PNI =0  DO 
540  PN“PNI 
J=0 

520  P  =  PN 

C  FIND  YO  VS  L  SLOPE  AND  CE  VS  L  SLOPE 
1  =  1 

IF(P.LT.UI)  )P-U  1  ) 

500  IF(P.LT.UIM))G0T0  510 
I  =  I>1 
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1  168 
1  169 
1  170 
117  1 
1  172 
1  173 
1  174 
1  175 
1  176 
1  177 
1  178 
1  179 
1  180 
1  181 
1  182 
1  183 
1  184 
1  185 
1  186 
1  187 
1  188 
1  189 
1  190 
1  19  1 
1  192 
1  193 
1  194 
1  195 
1  196 
1  197 
1  198 
1  199 
1200 
1201 
1202 

1203 

1204 

1205 

1206 

1207 

1208 
1209 
12  10 
12  11 
12  12 
12  13 
12  14 
12  15 
12  16 
12  17 
12  18 
12  19 
1220 
122  1 
1222 

1223 

1224 

1225 

1226 

1227 

1228 

1229 

1230 

1231 

1232 

1233 

1234 

1235 

1236 

1237 

1238 

1239 

1240 
124  1 


I F ( I  LT. ICT )G0T0  500 
P=L ( ICT) 

510  DIST*P-L(I) 

SSLP-6 . D0*CEE ( I .3) 

SLP=6 . D0*CEE ( I , 3 ) *DlST+2 . D0*CEE ( I , 2 ) 

I F ( DABS( SSLP-0 . DO ) .  LT  .  1 .D-10)G0T0  512 
C  THE  NEWTON-RAPHSON  METHOD 
PN=P-SLP/SSLP 
J  =  J+1 

I F ( J . LT . 100 ) GO TO  530 
IF(PNI  .LE.-4.D~2  )GOTO  550 
PN I =  PN 1-1 . D-2 
GOTO  540 

550  WR I TE ( 6 , 20 ) 

WR I TE ( 7 , 20 ) 

GOTO  560 

530  IF(DABS(PN-P) .GT. 1  D-5)G0T0  520 
C 

C  FIND  THE  TOTAL  AND  MAX.  COLLISION  EFFICIENCY. 

512  CET  =  ( Y0( ICT ) - Y0( 1 ) ) *DCOS( ALPHAR ) /THICK* 1 .D4 

CEMAX= ( (3 . DO* CEE ( I , 3 ) *DIST+2 . DO*CEE ( I , 2 ) ) *DIST+CEE ( I , 1 ) )*1 .D2 
. *DCOS ( ALPHAR ) 

WR I TE ( 6 , 10 )CEMAX . P , CET 
WR I  T  E  (  7  ,  1 0 )CEMAX , P , CET 
C 

C  DETERMINE  PLOTTING  PARAMETERS. 

560  IF ( LAYER. EO.  1  )CALL  PLTSZ(LS(  1 ) , LS ( I CT ) , YOS ( 1 ) , YOS ( I CT ) , 

. LPMIN , YOPMIN , PX , PY , I RX , I RY ) 

IF (LAYER. GT.  1  )CALL  PLTSZE(LS( 1 ) . LS ( I CT ) . YOS ( 1 ) , YOS ( I CT ) , 

. LPMIN, YOPMIN, PX , PY , IRX , IRY ) 

LS ( ICT+1 ) - LPM I N 

LSI  I  CT  +  2  )  =  XPAR ( 4 ,  I RX ) / 10.0*  *PX 
CALL  NEWPEN( 1 ) 

IF(Y0L  EO.O)GOTO  200 

C 

C  PLOT  THE  YO  VS  L  GRAPH. 

YOS (  ICT  +  1  )  =  YOPM I N 
Y  OS ( I CT  +  2 )  =  YPAR ( 4 , IRY ) / 10 . 0*  *  PY 
CALL  FACTOR! FACT ( YOL )*PLTFAC) 

CALL  ORIGIN! 999, 20.0. 13.0,5.0,5.0) 

CALL  AX2EP(XPAR(3. IRX), 3. 1+PX.O, 1.0) 

CALL  AX  I S2 (0.0. 0.0,  'L/C'  , -3. XPAR( 2, IRX ) , 0 . O , LPMI N , XPAR ( 4 , IRX )/ 
10.0*  *  PX , XPAR( 3 , IRX) ) 

CALL  A X  I S2 ( XPAR ( 2 . IRX ) ,0  0,  '  ' , -  1 . -YPAR (2 , I R Y ) , 90 . 0 ,  1  . 0 ,  1  .0, YPAR ( 3 

.  , IRY  )  ) 

CALL  AX2EP ( YPAR ( 3 ,  IRY  )  , 3, 1  +  PY .0,  1.1) 

CALL  AX  I S2( 0.0, 0.0,  ' YO/C '  ,4, YPAR(2, I RY ) , 90 . O , YOPM I N . YPAR ( 4 , IRY)/ 

.  10  0*  *PY , - YP AR ( 3,  IRY ) ) 

CALL  AX  I S2( 0.0. YPAR (2 .IRY),  '  '  ,  1 . -XPAR (2 , I RX ) . 0 . 0 ,  1  .  ,  1  .  . XPAR ( 3 . IRX 

.  )  ) 

C  PLOT  THE  YO  VS  L  POINTS 
CALL  LINEPfO. 15) 

CALL  LINE< LS. YOS. ICT .  1 . -  1 .0) 

F  *  1 

L  RG=LS ( ICT ) - L  S (  1 ) 

YORG  =  YOS (  ICT  )  - YOS (  1 ) 

DO  100  1=1,101 

LP (  I  )  =LS (  1  )  +  F  LOAT ( I  - 1 )/ 100 . O*  LRG 
120  IF( LP(  I ) . LE . LS( F+ 1  )  )GOTO  110 

F  =  F+  1 
GOTO  120 

110  D  =  LP( I  )-LS(F  ) 

Y  OP ( I  )  =  SNGL ( ( (CEE(F , 3)*D+CEE(F , 2  )  )  *D+CEE( F .  1 ) ) *D )+YOS< F ) 

100  CONTINUE 

Y0P(  102 )  =  Y0S(  ICT  +  1  ) 

Y0P( 103 ) = YOS ( ICT+2 ) 

LP(  102  )  =LS(  I C  T  +  1  ) 

L P (  103 )  =  LS(  ICT  +  2  ) 

C  PLOT  THE  YO  VS  L  LINE 

CALL  LINEUP.  YOP.  101 , 1.0.  1  ) 

C 

C  PLOT  THE  CF  VS  L  GRAPH 
200  I  F ( CE  L  EQ  0)GOTO  300 

CALL  FACTOR! FACT(CEL  )*PLTFAC  ) 
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1242 

1243 

1244 

1245 

1246 

1247 

1248 
1219 

1250 

1251 

1252 

1253 

1254 

1255 

1256 

1257 

1258 

1259 

1260 
1  26  1 
1262 

1263 

1264 

1265 
1  266 

1267 

1268 

1269 

1270 

127  1 

1272 

1273 

1274 

1275 

1276 
1  277 

1278 

1279 

1280 

128  1 
1282 
1  283 

1284 

1285 

1286 

1287 

1288 

1289 

1290 
129  1 
1  292 
1  293 

1294 

1295 

1296 
1  297 

1298 

1299 

1300 
1  30  1 
1  302 
1  303 
1  304 

1305 

1306 

13(W 

1  308 
1  309 
1  3  10 
13  11 
1312 
13  13 
13  14 
1315 


CALL  0RIGIN(999, 20 . 0 , 13.0,5.0.5.0) 

CALL  AX2EP ( XPAR ( 3 , IRX  )  ,3,  1+PX.O,  1.0) 

CALL  AX  I S2( 0.0. 0.0. 'L/C'  . -3 . XPAR<2 . IRX ) .0 . 0, LPMIN, XPAR( 4 . IRX )/ 

10 . 0*  *PX ,XPAR( 3. IRX)) 

CALL  AX  I S2 ( XPAR ( 2 . IRX ) .0  0.  '  ' , -  1 . -YPAR(2 , 10 ) . 90 . 0 . 0 . 0. 1  .0. YPAR(3 

.  10)  ) 

CALL  AX2EP(YPAR(3, 10). 3,0.0. 1 . 1 ) 

CALL  AXIS2(0  0.0  .0, 'COLLISION  EFFICIENCY  IN  % ' . 25 , YPAR ( 2 . 10 ) , 
.90.0.0.0. YPAR( 4, 10)* 10.0, -YPAR (3. 10)) 

CALL  AXIS2(0.0, YPAR ( 2 , 10) , '  ' , 1 , -20 . 0 . 0 . O . 1 . . 1 . .XPAR(3, IRX ) ) 

LRG  =  LS<  I CT ) - L  S (  1  ) 

F  =  1 

C  DETERMINE  PLOTTING  VALUES  OF  CE . 

DO  210  1=1,101 

LP( I  )  =  LS(  1 )  +  F  LOAT ( 1-1 )/lOO.O*LRG 
230  I F ( LP ( I ) . LE.LS(F+1 ) )GOTO  220 

F  =  F  +  1 
GOTO  230 

220  D  =  LP( I  )-LS(F  ) 

CEP ( I  )  =  SNGL ( ( 3  DO*CEE (F,3)*D+2. DO*CEE ( F , 2 ) )*D+CEE(F ,  1 ))* 100-0 
*COS<  SNGU  ALPHAR  )  ) 

210  CONTINUE 

LP(  102 )  =  LS( I CT  +  1  ) 

LP(  103 )  =  LS (  ICT+2  ) 

CEP ( 102 )=0.0 

CEP  (  103  )  =  YPAR ( 4 .  10 ) *  10 . O 
C  PLOT  THE  CE  VS  L  LINE . 

CALL  LINE1LP, CEP. 101, 1.0.1) 

300  IFfCEX. EO-O. OR. LAYER. GT. 1 )GOTO  400 
DO  310  KL  = 1 .NEL 
C 

C  PLOT  THE  CE  VS  X  GRAPH. 

IF( -LL(KL  )  . LE . L(  1  )  )GOTO  320 
310  CONTINUE 

320  DO  330  KU  = 1 ,NEU 

IF(LU(KU)  .GT.L( ICT))GOTO  340 
330  CONTINUE 

340  XRG  =  SNGL( XL ( KL  )  +  XU( KU  )  ) 

XLF  =SNGL \ -XU(KU) ) 

C0  =  0 
I  I  =ICT-  1 

DO  350  KK= 1 , 201 
X=XLF+XRG/200 .  * F  LOAT ( KK  -  1  ) 

XP ( KK  )  =X 

C  DETERMINE  VALUE  OF  L  FOR  EACH  X. 

IF( X .GT .0. )GOTO  360 

CALL  SFC(DBLE( -X) . Y , 1 . 1 , ZZ) 

GOTO  370 

360  CALL  SFC(DBLE(X),Y.O. 1 ,ZZ) 

ZZ=-ZZ 

370  I F ( CO . EQ  1  ) GO T 0  380 

I F ( ZZ . GT . L (  ICT)  )GOTO  380 
IE ( ZZ  GT . L(  I  I) )GOTO  4 10 
11=11-1 

I F (  I  I  . EQ.O)GOTO  390 
GOTO  370 
390  C0=1 

380  CE P ( KK  )  =0  O 

GOTO  350 

4  10  D  =  7  Z  -  L (  I  I  ) 

CEP( KK ) - SNGL ( ( 3  00*CEE( 1 1 ,3)*0  +  2 .DO*CEE( 1 1 . 2 ) >*D+CEEU I ,  1 ) ) *  100 
*  COS ( SNGL ( ALPHAR  )  ) 

350  CONTINUE 

C  DETERMINE  THE  PLOTTING  PARAMETERS. 

CALL  PLTSZE ( XP<  1 ) . XP <  20 1 ) . O . DO . 99 . 9 , XPMIN , CEPMIN . PX , P Y . I RX . I RY ) 
XP( 202 )  =  X PM  I N 

XP( 203 ) -XPAR1 4 . IRX )/ 10  O*  *PX 
CEPf  202 1=0.0 

C E P ( 203  )  =  YPAR ( 4 .  10)*  10  0 

r  pi  or  cf  vs  x  axes 

CALL  FACT OR  1  FACT (CEX  )*PLTFAC  ) 

CALL  OR  I G I N( 999 . 20  0.  13.0.5  0.5  O) 

CALL  AX2EPI  xpAR( 3.  IRX  )  .3.  1+PX.O.  1  O) 

CALI  A  X  I S  7 ( O  0.0  O.  ' X/C '  .  3 . XPAR( 2 .  IRX ) .0  0. XPMIN. XPAR( 4 ,  IRX ) 


1316  . / 10 .0* *PX , XPAR( 3 . IRX ) ) 

1317  CALL  AX  I S2 ( XP AR ( 2 . IRX) ,0.0. '  ' , -1.-YPARI2. 10) ,90  0,0.0, 1  0.YPARO, 

1318  .  10)  ) 

1319  CALL  AX2EP(YPAR(3. 10) .3.0.0. 1 . 1 ) 

1320  CALL  AXIS2(0. 0,0.0, 'COLLISION  EFFICIENCY  IN  % '  , 25 , YPAR ( 2 .  10 ) . 

132  1  .90.0,0.0, YPAR ( 4 .  10) *  10. O. - YPAR ( 3 ,  10)  ) 

1322  CALL  AXIS2(0.0,YPAR(2. 10), '  1 , -20  0,0.0, 1 ., 1 . ,XPAR(3. IRX) ) 

1323  C  PLOT  THE  CE  VS  X  LINE. 

1324  CALL  LINE( XP.CEP.201 . 1.0. 1) 

1325  400  RETURN 

1326  END 

1327  C 

1328  C 

1329  SUBROUTINE  PLTSZ( XMI N , XMAX , YM I N , YMAX . XL . YB , PX . PY .  I RX , I RY  ) 

1330  C 

1331  C  WRITTEN  BY:  M  OLESKIW  ON  800627  LAST  MOD  I F I E D : 80 1 0 1 8 

1332  C 

1333  C  DETERMINE  PARAMETERS  NECESSARY  FOR  SCALING  OF  A  PLOT  AND  ITS  AXES 

1334  C 

1335  REAL  XPAR(4,  10),YPAR(4,  1 0 ) . XD . F LOA T , A  I  NT , XM I N . XMAX . 

1336  .XL. YD, YMIN, YMAX , YB , DX ,DY , XR , YT 

1337  C 

13  38  INTEGER  PX , PNX ,PY , PNY .  I . J .  I  X . I RX .  I  NT ,  I Y .  I RY .  I F I  X 

1339  C 

1340  COMMON/ PL TP RM/ XPAR.YPAR 

134  1  C 

1342  C  IN  XM I N= 

1343  C  IN  XMAX= 

1344  C  IN  YM IN= 

1345  C  IN  YMAX  = 

1346  C  OUT  XL=LEFT  EDGE  OF  PLOT 

1347  C  OUT  YB  =  BOT T OM  EDGE  OF  PLOT 

1348  C  OUT  PX=POWER  OF  TEN  IN  X-AXIS  RANGE 

1349  C  OUT  P Y  =  POWER  OF  TEN  IN  Y-AXIS  RANGE 

1350  C  OUT  I RX  =M I N .  LENGTH  OF  X  AXIS- 

1351  C  OUT  I R Y  =M I N .  LENGTH  OF  Y  AXIS. 

1352  C 

1353  10  F0RMAT18F 10.0) 

1354  C 

1355  C  READ  IN  PLOTTING  PARAMETERS 

1 356  DO  101  I  - 2 .  10 

1357  RE  AD (3.  10 )( XPAR( J. I  ) ,0-1 .4) , (YPAR(J, I  )  . U- 1  . 4  ) 

1358  101  CONTINUE 

1 359  C 

1360  ENTRY  PLTSZE ( XMIN , XMAX . YMIN . YMAX , X!  . YB . PX , PY ,  I RX .  I RY  ) 

1361  PNX  =0 

1362  PNY  =0 

1363  C 

1364  C  DETERMINE  THE  PLOTTING  RANGE  OF  THE  X  VARIABLE 

1365  100  PX -PNX 

1366  XD= ( XMAX- XMIN )♦ 10 . 0* *PX 

1367  IF(XDGT.  10.0 ) PNX  =  PNX -  1 

1 368  IF(XD  LT.  1  00001  )PNX  =  PNX+ 1 

1369  I F ( PNX . NE  PX  )GOFO  100 

1370  C  PX  GIVES  1 / ( POWER  OF  TEN)  OF  THE  X  VARIABLE  PLOTTING  RANGE 

1371  I  X  =  1 

1372  120  I R  X  = I  NT ( XD  )  +  I  X 

137  3  DX  =  FLOAT ( IRX ) / 10 . 0* *PX/XPAR(  1 ,  I RX  ) 

1374  C  SET  THE  X  VALUE  AT  THE  LEFT  GRAPH  EDGE 

1375  IF( XMIN. LT  0)XL=AINT< XMIN/DX- 1  O) *DX 

137  6  I F<  XMIN . GE  . 0)XL  =  AINT< XMIN/DX ) *DX 

1377  XR~XL*XPAR(  1  ,  IRX  )*DX 

1378  IF ( XR  GE . XMAX  AND  IRX  NE  3  AND 

1379  IRX  NE  6  AND  I RX  NE  7  AND  I RX  NE  9  ) GOTO  105 

1 380  I  X  =  I  X  4  1 

1381  GOTO  120 

1382  105  I F (  I F  I  X (  ( XR  xMAX  1/DX)  LE  I  F  I  X  C  (XMIN-XL )/DX )  ) GO  T  O  110 

1383  C  CENTRE  THE  PLOT 

1384  XL-XL-DX 

1385  XR=XR-DX 

1386  GOTO  105 

1 387  C 

1388  C  DETERMINE  THE  PLOTTING  RANGE  OF  THE  Y  VARIABLE 

1 389  1  10  P Y  =  PNY 


1390 

1391 

1392 

1393 

1394 

1395 

1396 

1397 

1398 

1399 

1400 

1401 

1402 

1403 

1404 

1405 

1406 
1  407 

1408 

1409 
14  10 
14  11 
14  12 
14  13 
14  14 
14  15 
14  16 
14  17 
14  18 
14  19 
1420 

142  1 

1422 

1423 

1424 

1425 

1426 

1427 

1428 

1429 

1430 

143  1 

1432 

1433 

1434 

1435 

1436 

1437 

1438 

1439 

1440 

144  1 
1442 
144  3 

1444 

1445 

1446 

1447 

1448 

1449 

1450 
14  5  1 

1452 

1453 
1  454 

1455 

1456 
1  457 
1458 
1  459 
1460 
1  46  1 

1462 

1463 


YD=  ( YMAX-YMIN)* 10 . 0*  *PY 
IF( YD.GT .9 . 99999 )PNY=PNY- 1 
I F ( YD  LT .  1  . 0 ) PNY  =PNY  + 1 
I F ( PNY  NE  . PY  )G0T0  110 

C  PY  GIVES  1 / ! POWER  OF  TEN)  OF  THE  Y  VARIABLE  PLOTTING  RANGE 
I  Y=  1 

130  IRY  =  INT< YD  )  +  IY 

DY  =  F  LOA  T (  IRY ) / 10  O* * PY/YPAR (  1  , IRY  ) 

C  SET  THE  Y  VALUE  AT  THE  BOTTOM  OF  THE  GRAPH 
I F ( YMIN  LT  0  0)YB=AINT( YMIN/OY-1 .0)*DY 
IFfYMIN  GE  00)YB=AINT( YMIN/DY ) *DY 
Y  T  =  v'B  +  v  PAR  (  1  .  I  RY  )  *DY 
I F ( YT  GE  YMAX  )GOTO  135 
I  Y  =  2 


GOTO  130 

135  I F 1  I F  I  X (  ( V  T- YMAX )/DY )  . LE  I F I  X ( ( YMIN-YB )/DY ) )GOTO  140 
C  CENTRE  THE  PLOT 
YB= YB -DY 
YT=YT -DY 
GOTO  135 
140  RETURN 
END 


C 

C 


SUBROUTINE  ICING! CETOL , I CE . BOTH . F A  I L  ) 


C 


C  WRITTEN  BY  M  OLESKIW  ON : 8007 1 3  LAST  MODI F I  ED : 80 1 227 


C 

C  CALCULATE  AMOUNT  0 r  ACCRETION  AND  DETERMINE  A  NEW  SET  OF  AEROFOIL 
C  SURFACE  ELEMENT  ENDPOINTS  AFTER  DETERMINING  THE  AEROFOIL 
C  NOSE  LOCATION 

C 

DOUBLE  PRECISION  XN, VN, XNN , YNN , XUR( 101  )  , YUR( 101 ) , 

CU( 100. 3) .CL ( 100. 3) . XLR( 101 ) , YLR( 10 1 ) , L ( 5 1 ) , Y0( 5 1 ) . 

D. CEE! 50. 3 ), CETOL . K . DSIGN, XLRN( 101 ) . YLRN( 101 ) . 

S30.C30.NSURF , XURN(  101 ) , YURN(  101  )  . CEU(  101).CEL(101),D1.D2, 

XUT (  10  1  )  . XLT (  101  )  . XU(  101  )  , YU(  101  )  . XL(  101 ) . YL(  101 ) . 

YUT{  101  )  .  YLT(  101 ) .DABS . DSQRT . LU( 101 ) . LL< 101 ) , ICE . 

PP . TOL . LE , RE , I  CEE . ALPHAR , DCOS 

DOUBLE  PRECISION  NSURF Y , XRMI N . XNP , YNP , XURTLP . XLRTLP 

INTEGER  BOTH, J . NCOU , NCOL . NCOUN . NCOLN , I CT , I CU . I CL . I , I ER . NOAC . ONCE , 

IM,  I  US.  ILS ,  IK . FAIL .RUN, NEU ,NEL . I U ( 5 1  ) .  I L ( 5 1 ) , I UN( 5 1 ) , I LN( 5 1 ) , 

I  1  .  12 . J  1  . J2 .KK.KL . LLL . IXU<  101 ) . IXL(  101 ) . IZU(  101 ) . IZL(  101 ) , IUU, ILL 

C 

COMMON  ALPHAR /AER03/NC0U . NCOL /NOSE / XN . YN/ FO I L / XUR , YUR , 

. XLR . YLR/R0TP/C30. S30/CEM/PP/ I ND/NSURF Y . I CEE . I , J . RUN/ AER04/NEU . 
NEL/COL/L . YO, ICT, ICU, ICL/EFF/CE E/SFCS/XU . YU . XL . YL/LG/LU , LL 
/SPLINE/CU. CL /ENDS/ I U . I L /NNOSE / XNP . YNP , XURTLP , XLRTLP 

C 

EXTERNAL  NSURF 


C 

C  IN  CETOL-CRITERI ON  FOR  DETERMINING  THE 
C  SEGMENT  ENDPOINTS 

C  IN  ICE-MAX  THTCKNESS  OF  ICE  ACCRETION 
C  IN  BOTH= TRAJECTORIES  FOR  BOTH  SFCS  (O 
C  OUT  FAIL=FAI LURE  INDICATOR 
C 

10  FORMAT! '  FAILURE  TO  CONVERGE  TO  NEW 
20  FORMAT ( ' 1 ENDPT  X  COORD  Y  COORD 
30  FORMAT! '  ' . F 14 . 5 . F 10 . 5 . F 1 7  5.F12  4) 

40  FORMAT!  '  '  ) 


NEED  FOR  NEW  CONTROL 

! ASSUMING  CE= 100%) . 
OR  1  ) 


NOSE  POSITION' ) 

DIST  FROM  NOSE  COLL . 


EFF 


) 


C 

XURTLP=XUR! NEU ) 

XLRTLP  =  XLR( NEL  ) 

J  -  I  CL 
NOAC  =0 
ONC  E  =0 
C 

C  FOR  THE  UPPER  SFC 

DO  100  1=1. NFU 
I F ! NOAC . EQ  1 IGOTO  1 15 

C  DETERMINE  THE  APPROPRIATE  CE  VS  L  SEGMENT 
110  IF(L.U(I)  lE.UJMMGOTO  120 

J  =  J  +  1 


50 


1464 

1465 

1466 

1467 

1468 

1469 

1470 
147  1 

1472 

1473 

1474 

1475 

1476 

1477 

1478 

1479 

1480 


1487 

1488 

1489 

1490 

1491 

1492 

1493 

1494 

1495 

1496 

1497 

1498 

1499 

1500 

1501 

1502 

1503 

1504 

1505 

1506 

1507 

1508 

1509 
15  10 
15  11 
15  12 
15  13 
15  14 
15  15 
15  16 
15  17 
15  18 

1519 

1520 

152  1 

1522 

1523 

1524 

1525 

1526 

1527 

1528 

1529 

1530 

153  1 

1532 

1533 

1534 

1535 

1536 

1537 


I F ( J . L  T . ICT) GO  T  0  110 
N0AC=1 

C  NO  ACCRETION  REGION  ON  TOP  SFC . 

115  CEU( I  )  =0 . DO 

XURN( I )  =XUR ( I ) 

YURN( I ) = YUR ( I ) 

GOTO  100 

U°  CEUl"  I  )  =  (  (  3^D0*CEE  ( d ,  3 )  *0  +  2  .  DO*CEE (  d,2)  )*D+CEE(d,  1  )  )  +  DCOS  (  ALPHAR  ) 

I  F  ( DABS ( CU( I ,  1  )  )  .  LT  .  1  .D-20)G0T0  150 
K  =  -  1  .00/CU( r .  1 ) 


XURN ( I ) =XUR ( I )+DSI GN( DSQRT ( ICE*ICE*CEU( I ) *CEU (  I )/(  1  DO+K*K ) ) , K ) 
YURN( I ) =  YUR ( I )+K*  ( XURN( I )-XUR( I ) ) 

GOTO  100 

C  GROWTH  IN  Y  AXIS  DIRECTION 


1481 

150 

XURN< I )  =  XUR ( I ) 

1482 

YURN( I )  =  YUR ( I ) +CEU( I ) *  ICE 

1483 

100 

CONTINUE 

1484 

DO  160  1=1 .NCOU 

1485 

I UN( I )=IU( I) 

1486 

160 

CONTINUE 

NCOUN=NCOU 


CHECK 


335 


C  CHECK 

C  CHECK 
325 


320 

330 


340 

350 


C  CHECK 
3  15 


360 

370 

C 

C  SHIFT 


380 


300 

390 


180 


FOR  NEED  OF  CREATING  NEW  CONTROL  ENDPTS.  ON  UPPER  SFC 

DO  300  1=2, NCOU 

I F ( ONCE . EO.O)GOTO  335 

ONCE  =0 

GOTO  300 

I  1  =  I UN( I ) 

12=1 UN (1-1 ) 

I F ( I  1 .EQ. 12+1 )GOTO  300 

IF(CEU(I2)  .EQ.O. DO ) GO TO  390 

FOR  ZERO  CE  BETWEEN  CONTROL  ENDPTS. 

IF(CEU( I  1 ) .EQ.O. DO ) GOTO  330 
FOR  RAPID  CHANGE  IN  CE 

I F ( DABS ( CEU ( I  1 ) -CEU( I  2 ) ) . LT . CETOL )GOTO  315 
J 1  =  I  2+  1 
J2=  I  1 

DO  320  d  =  d1 , J2 

IF(DABS(CEU(d)-CEU( J-1 ) ) . GE .CETOL/ 1 .2DO)GOTO  350 
CONTINUE 
GOTO  360 
d  1  =  I  2+  1 
d2  =  I  1 

DO  340  d=d 1.02 

IF(CEU(d) .EQ.O. DO .AND.CEU(d~1 ) .GE .CETOL/2. DO )GOTO  350 
CONTINUE 
GOTO  325 
KK=d' 1 

IF(d.EQ.d1 ) KK - d 
GOTO  370 

IF  DISTANCE  BETWEEN  CONTROL  ENDPTS-  IS  INCREASING  SUBSTANTIALLY 
D1=DS0RT( (XUR( I  1 ) -XUR( 12) ) *  *2+ ( YUR ( I  1 ) - YUR ( I  2 ) ) *  * 2 ) 

02  =DSQRT ( ( XURN ( I  1 )-XURN( 12) ) *  * 2+ ( YURN( I  1 )-YURN(  12)  )**2) 

I F( D2 . LT .  1  . 25D0*D1  )GOTO  300 
KK= ( I  1  +  12  )/2 
ONCE  = 1 

KL  =NCOUN- I +  1 


INDICES  OF  CONTROL  ENDPTS. 

DO  380  LLL  = 1 .KL 
IUN( NC0UN+2-LLL ) = IUN( NCOUN+ 1 
CONTINUE 
NCOUN=NCOUN+ 1 
I  UN (  I  )=KK 
CONTINUE 
=  1 

DO  170  1  =  1  , NEU 
I F ( IUN(d)  .EQ.  I  )GOT  0  180 
IXU<  I  )=0 
GOTO  170 
I XU( I  ) *  1 


TO  MAKE  ROOM  FOR  A  NEW  ONE 


-LLL) 


51 


r— r 


i 


s 


► 

K 

h . 


1538 

1539 

1540 

1541 

1542 

1543 

1544 

1545 

1546 

1547 

1548 

1549 

1550 

1551 

1552 

1553 

1554 

1555 

1556 

1557 

1558 

1559 

1560 

156  1 

1562 

1563 

1564 

1565 

1566 

1567 

1568 

1 569 

1570 

157  1 


170 


190 

195 


J  =  J+  1 
CONTINUE 
WR I TE ( 7 , 20 ) 

DO  190  1*1, NEU 

WRITE ( 7.30)XU( I ) , YU( I ) ,LU( I ) ,CEU( I ) 
I  F ( CEU ( I  )  . EQ . O . DO )GOTO  195 
CONTINUE 

I F ( BOTH  EQ . O )GOTO  590 
J=ICL+  1 
NO AC  =  0 


FOR 


VS  L 
220 


SEGMENT . 


THE  LOWER  SFC . : 

DO  200  1*1 .NEL 
I F ( NOAC . EQ .  1  )GOTO  215 
C  DETERMINE  THE  APPROPRIATE  CE 
2  10  I  F ( ~LL( I  )  .GT . L( J) )GOTO 

J  =  J-  1 

I F ( J . GT . 0 ) GOTO  210 
NOAC  *  1 

C  NO  ACCRETION  REGION  ON  LOWER  SFC. 

215  CEL ( I ) =0 . DO 

XLRN( I )=XLR( I ) 

YLRN( I  )  = YLR (  I  ) 

GOTO  200 
D*-LL(I)-L(J) 

CEU I  )  =  ( ( 3  D0*CEE ( J  t  3 ) *D  +  2  DO*CEE ( J , 2 ) ) *D+CEE  (  J  .  1  )  )  *DCOS ( ALPHAR ) 
I F ( DABS ( CL ( I .  1  )  )  . LT .  1 . D-20)G0T0  250 
K*-1 . DO/CL ( I , 1 ) 


220 


NEW  ENDPTS . : 

XLRN(  I  )  =  XLR  (  I ) -DS IGN( DSORT ( ICE  MCE *CEL ( I  )  *CEL (  I  )/(  1  D0+K*K ) ) .K) 
YLRN( I )  =  YLR ( I )+K*(XLRN( I )-XLR( I ) ) 

GOTO  200 

GROWTH  IN  Y  AXIS  DIRECTION 


1572 

250 

XLRN( I ) =XLR (  I  ) 

1573 

YLRN( I  ) =  Y L R ( I  ) -CEL (  I )*ICE 

1574 

200 

CONTINUE 

1575 

DO  260  1=1, NCOL 

1576 

I LN( I  )  =  I  L (  I  ) 

1577 

260 

CONTINUE 

1578 

NCOLN  =  NCOL 

1579 

ONCE  =0 

1580 

C 

158  1 

C  CHECK  FOR  NEED  OF  CREATING  NEW  CONTROL  ENDPTS.  ON  LOWER  SFC. 

1582 

DO  400  1=2, NCOL 

1583 

I F ( ONC E . EQ  0 )GOTO  435 

1584 

ONCE  =  0 

1585 

GOTO  400 

1586 

435 

I  1  =  ILN(  I  ) 

1587 

I  2  =  I LN (  H  ) 

1588 

IF(11  EQ. 12+1 )GOTO  400 

1589 

IF(CEL( 12) . EQ . 0 . DO ) GOTO  905 

1590 

C  CHECK  FOR  ZERO  CE  BETWEEN  CONTROL  ENDPTS. 

1591 

I F ( CEL ( I 1 1 .EQ.O. DO ) GO TO  430 

1592 

C  CHECK  FOR  RAPID  CHANGE  IN  CE . 

1593 

425 

IF (DABS ( CEL ( I1)-CEL(I2)).LT  CETOL )GOTO  415 

1594 

U 1  *  1 2+1 

1595 

J2  =  I  1 

1596 

DO  420  U  =  J1 . J2 

1597 

IF(DABS(CEL(U)-CEL( J-1 ) )  GE  .CETOL/ 1 .200) GOTO  450 

1598 

420 

CONTINUE 

1599 

GOTO  460 

1600 

430 

J  1  =  I  2  +  1, 

1601 

U2=I  1 

1602 

DO  440  J=J1 . J2 

1603 

IF (CEU  J  ). EQ.O. DO. AND . CEL( J- 1 )  GE . CETOL/ 2 . DO) GOTO  450 

1604 

440 

CONTINUE 

1605 

GOTO  425 

1606 

450 

KK  =  J -  1 

1607 

IMU.EQ.J1  )KK  =  J 

1608 

GOTO  470 

1609 

C  CHECK  IF  DISTANCE  BETWEEN  CONTROL  ENDPTS  IS  INCREASING  SUBSTANT: 

16  10 

4  15 

0 1 =DSQRT ( <  XLR( I  1  )-XLR< 1 2  )  ) * *2*< VLR( 1 1 ) - YLR(1 2 ) ) *  *  2 ) 

16  1  1 

02  =  DSORT ( ( XLRN( I 1 ) -XLRN( I 2 ) ) ♦ ♦ 2  + ( YLRN( I 1 )-YLRN( 12) )**2 ) 

52 


16  12 

IF(D2.LT  1  25DO*D1)GOTO  400 

1613 

460 

KK * ( I  1  +  I 2  )  /2 

1614 

ONC  E  = 1 

16  15 

4  70 

KL -NCOLN- I + 1 

16  16 

C 

16  17 

C  SHIFT  INDICES  OF  CONTROL  ENDPTS.  TO  MAKE  ROOM 

FOR 

A  NEW  ONE . 

16  18 

DO  480  LLL= 1 , KL 

16  19 

I LN( NCOLN  +  2 - LLL )  =  I LN( NCOLN+ 1  - LLL ) 

1620 

480 

CONTINUE 

162  1 

NCOLN =NCOLN+ 1 

1622 

I LN ( I ) =KK 

1623 

400 

CONTINUE 

1624 

905 

J=1 

1625 

DO  270  1=1 .NEL 

1626 

I F ( I LN (  »)  .EQ. I )GOTO  280 

1627 

I  XL ( I  )  =0 

1628 

GOTO  270 

1629 

280 

I  XL (  I  )  =  1 

1630 

J  =  J+1 

163  1 

270 

CONTINUE 

1632 

WR I T  E ( 7 , 40 ) 

1633 

DO  230  1=1 .NEL 

1634 

WRITE ( 7 , 30 )XL( I  ),YL(l),LL(I  )  ,CEL( I ) 

1635 

I F ( CE L ( I  )  EQ.O. DO ) GOTO  900 

1636 

230 

CONTINUE 

1637 

GOTO  900 

1638 

C 

1639 

C  UPPER  &  LOWER  SFCS.  MIRROR  IMAGES;  NOSE  STAYS 

ON 

THE  X-AXIS. 

1640 

590 

DO  595  1  =  1  , NEU 

164  1 

XLRN( I )=XURN( I ) 

1642 

Y LRN( I  )  =  -YURN( I ) 

1643 

I  XL (  I )  =  IXU(  I  ) 

1644 

595 

CONTINUE 

1645 

GOTO  930 

1646 

C 

1647 

C  FIND  NEW  NOSE  LOCATION  USING  THE  GOLDEN  SECTION  SEARCH  METHOD 

1648 

C  OF  DETERMINING  THE  MIN.  VALUE  OE  THE  NEW  SURFACE  X-COORD. 

1649 

900 

I  CE  E  =  I CE 

1650 

RUN  =  0 

1651 

J  =  1 

1652 

1  =  1 

1653 

DO  9 10  KK  = 1  , NCOL 

1654 

IF(LL(KK) .  GE  . -PP )GOTO  920 

1655 

9  10 

CONTINUE 

1656 

920 

TOL  =  1  . D - 5 

1657 

F  A  I L  =0 

1658 

LE  =  1  . D-  10 

1659 

RE  =  XLR ( KK  ) 

1660 

CALL  ZXGSN(NSURF,LE,RE,T0L,XRMIN. IER) 

1661 

IF(  IER . LT .  129 .OR . IER -GT .  132  )GOTO  950 

1662 

F  A  I  L  =  1 

1663 

WR I TE ( 6 , 10) 

1664 

WR I T  E ( 7 ,  10) 

1665 

GOTO  720 

1666 

C  NEW 

NOSE  COORDS . : 

1667 

950 

YNN=NSURFY 

1668 

XNN=NSURF( XRMIN) 

1669 

C 

1670 

C  DE- 

ROTATE  NEW  UPPER  8.  LOWER  SFCS-  ABOUT  PREVIOUS 

NOSE  POSITION 

167  1 

930 

DO  500  1=1 ,NEU 

1672 

XUT(  I  )  =  XURN(  I  )*C30YURN( I  )*S30+XN 

1673 

YUT(  I  )=XURN(  I  ) ♦ S30+  YURN(  I  )*C30+YN 

1674 

500 

CONTINUE 

1675 

DO  510  1=1 .NEL 

1676 

XL  T ( I ) - XLRN ( I  ) 'C30+YLRNI I )*S30+XN 

1677 

YLT( I  )  =  - X LRN( I  )*S30+YLRN( I  )*C30+YN 

1678 

5  10 

CONTINUE 

1679 

I F ( BOTH . EO . 1 )GOTO  520 

1680 

XNN=XUT( 1 ) 

168  1 

YNN=YUT( 1 ) 

1682 

I  M=  1 

1683 

520 

XU( 1 ) =XNN 

1684 

X  L (  1  )  =  XNN 

1685 

YU(  1  )  =  YNN 

1686 

YL( 1 ) = YNN 

1687 

IUU  =  1 

1688 

I  LL=  1 

1689 

IF(BOTH. EO.O)GOTO  625 

1690 

C 

169  1 

C  SEE 

IF  ANY  LOWER  SFC.  ENDPTS.  ARE  ABOVE  THE  NEW 

NOSE 

POSITION 

1692 

C  & 

THUS  BELONG  ON  THE  UPPER  SFC. 

1693 

DO  610  IM  = 1 , NEL 

1694 

IF (DABS ( YLT(IM)-YNN) . LT. 1 ,D-4)G0T0  620 

1695 

IF(YLT(IM).LT. YNN) GOTO  630 

1696 

6  10 

CONTINUE 

1697 

620 

I F ( IM . GT . 2 )GOTO  640 

1698 

I F ( I M . EQ . 2 )GOTO  650 

1699 

C  SAME  NOSE  INDEX 

1700 

625 

I  US  =  2 

1  701 

I  LS  =  2 

1702 

GOTO  665 

1703 

C  NEW 

NOSE  IS  NEAR  FIRST  ENDPT .  BELOW  PREVIOUS  NOSE 

1704 

650 

IUS=  1 

1705 

I  L  S  =  3 

1706 

GOTO  665 

1707 

C  NEW 

NOSE  IS  NEAR  SECOND  OR  GREATER  ENDPT.  BELOW 

PREVIOUS  NOSE 

1708 

640 

I K  =  I M  -  2 

1709 

DO  670  1=1 . IK 

1710 

IUU=IUU+ 1 

17  11 

XU( I UU ) =XLT ( IM- I ) 

17  12 

YU ( IUU  )  =  YLT (  IM- I ) 

17  13 

IZU(  IUU)  =  IX‘_(  IM-I  ) 

17  14 

670 

CONTINUE 

17  15 

IUS=  1 

1  7  16 

1LSMM+1 

17  17 

665 

IZU(  1  )=  1 

1  7  18 

IZL(  1  )=1 

1  7  19 

GOTO  660 

1720 

630 

I F ( I M . GT . 2 ) GOTO  680 

172  1 

C  NEW 

NOSE  IS  BETWEEN  FIRST  &  SECOND  ENDPTS.  ON  LOWER 

SFC  . 

1722 

IUS=  1 

1723 

I  LS  =  2 

1724 

GOTO  666 

1725 

C  NEW 

NOSE  IS  BELOW  SECOND  ENDPT.  ON  LOWER  SFC. 

1726 

680 

I  K  =  I M  -  2 

1727 

DO  690  1  =  1  .  IK 

1728 

IUU=IUU+1 

1729 

XU(  I UU )  =  X  L  T ( IM-I  ) 

1730 

YU ( IUU)=YLT( IM-I ) 

173  1 

IZU( I UU  )  =  I  XL (  IM-I ) 

1732 

690 

CONTINUE 

1733 

IUS=1 

1734 

I LS  =  I M 

1735 

666 

I  ZU (  1  )  =  1 

1736 

I ZL (  1  )=1 

1  737 

660 

DO  700  I  =  I  US , NEU 

1738 

I UU= I UU+ 1 

1739 

XU( I UU ) =XUT ( I ) 

1  740 

YU (  IUU  )  =  YUT( I  ) 

174  1 

IZU( IUU)  =  IXU<  I  ) 

1742 

I F (  I  E  0  I  US  AND  IUU . LT . 3  )  I ZU( IUU ) =0 

1  743 

700 

CONTINUE 

1744 

DO  710  I  =  I L  S , NE  L 

1745 

I LL  =  I LL  +  1 

1  746 

X L ( ILL  )=XLT(  I  ) 

1747 

YL(  ILL  )  = YLT (  I  ) 

1748 

I ZL ( I LL  )  =  I XL ( I  ) 

1  749 

7  10 

CONTINUE 

1750 

NEU= IUU 

175  1 

NEL  =  I  LL 

1752 

XNP=XN 

1753 

YNP  =  YN 

1  754 

XN  =  XNN 

1  755 

YN=  YNN 

1756 

IUU-  1 

1  757 

DO  730  1=1. NEU 

1  758 

I  F (  IZU<  I  > . EO . 0 )GOTO  730 

1  759 

IU( IUU)=I 

760 

761 

762 

763 

764 

765 

766 

767 

768 

769 

770 

771 

772 

773 

774 

775 

776 

777 

778 

779 

780 

78  1 

782 

783 

784 

785 

786 

787 

788 

789 

790 

79  1 

792 

793 

794 

795 

796 

797 

798 

799 

800 
801 
802 

803 

804 

805 

806 

807 

808 
809 
8  10 
8  1  1 
8  12 
8  1  3 
8  1  4 
8  15 
8  16 
8  1  7 
8  18 
8  19 
820 
82  1 
822 

823 

824 

825 

826 

827 

828 
829 
8  30 
83  1 

832 

833 


IUU=IUU+1 
730  CONTINUE 

I  LL  =  1 

DO  740  1=1 , NEL 
I F ( I ZL ( I ) . EQ . O )GQTO  740 
I  L ( I LL  )  =  I 
I Ll  =  I LL+  1 
740  CONTINUE 

NCOU= I UU- 1 
NCOL  = I LL- 1 
720  RETURN 
END 
C 

c 

SUBROUTINE  GROWTH( I CEPLA , LYRMAX , PLTF AC . TRJPL A ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON : 8007 1 3  LAST  MOD I F I  ED : 80 1022 

C 

C  PLOTS  SUCCESSIVE  AEROFOIL  OUTLINES  WITHIN  VIEW  WINDOW 

C 

REAL  XGR ( 104, 10).YGR( 104, 10 ) , PLTF AC , XMIN , XMAX , YMI N . YMAX , 

. XPLT ( 104) , YPLT( 104 ) ,XGRE( 103, 10) , YGRE( 103, 10) .XPLTE( 101 ) , 

. YPLTE (  101  ) 

C 

INTEGER  I T ( 10) , XZ , Y Z . LYRMAX , ICE PL A, ITT. I , J , TR JPL A , L YRM 1 , 

. I TE ( 10) , ITTE 
C 

COMMON/GROW/XGR, YGR.XGRE . YGRE . ITE . IT/GRID/XMIN , XMAX , YMIN, YMAX . XZ , 
.  YZ 
C 

C  IN  I CE PL A -PLOT  ACCRETION  OUTLINE.  (O  OR  1) 

C  IN  LYRMAX=NO .  OF  LAYERS  TO  BE  ACCRETED. 

C  IN  PLTFAC=PLOT  EXPANSION/REDUCTION  FACTOR. 

C  IN  TR JPLA=PLOT  TRAJECTORIES.  (O  OR  1) 

C 

I F ( ICEPLA.E02 )GOTO  120 
C  DRAW  AXES. 

CALL  NEWPEN( 1 ) 

CALL  ORIGIN( 999 , 2 1 .0,10.5,5.0,5.0) 

CALL  AX2EP(3.5.3,2,0,09) 

CALL  AX  I S2 ( 0 .  .0  ,  ' X/C  '  .-3,21.  ,0., XMIN, ( XMAX -XMI N ) /2 1  ,3.5) 

CALL  AXIS2(21  .  ,0.  ,  '  ' . - 1 . -  10 . 5 . 90 . .0 . .0 . . 1 . 75 ) 

CALL  AX2 E P (  1 .75,3,3,0.  1.1) 

CALL  AXIS2(0.  ,0.  ,  ' Y/C'  ,3,  10.5,90.  ,YMIN.  ( YMAX - YM I N ) / 1 0 . 5 , -1.75) 
CALL  AXIS2(0.  ,  10.5,  '  '  .  1 , -2  1  .  , 0  .  , XMIN ,  ( XMAX -XMI N ) /2 1  .,3.5) 

L YRM 1 =LYRMAX+ 1 
120  DO  100  1=1, L YRM 1 

I T T= I T ( I  ) 

I TTE  =  I TE ( I  ) 

DO  1 10  J= 1 . ITT 
XPLT ( J ) =XGR ( J . I ) 

YPLT ( J ) =YGR ( J , I ) 

110  CONTINUE 

DO  2  10  J=  1  , ITTE 
XPLTE(J) =XGRE ( J , I ) 

YPLTE(J) = YGRE ( J , I ) 

210  CONTINUE 

CALL  NE WPEN( 3 ) 

C  DRAW  ACCRETION  OUTLINES. 

CALL  LINE (XPLT, YPLT ,IT(I)-2, 1.0,0) 

CALL  L I NEP ( O . 1 ) 

C  PLOT  CONTROL  SEGMENT  ENDPTS . 

CALL  LINE(XPLTE, YPLTE, ITE( I )~2. 1,-1 .0) 

100  CONTINUE 

RETURN 
END 
C 
C 

SUBROUTINE  TRAJEC( TYPE . TRUPLA . TH I CK . A T . BOTH ) 

C 

C  WRITTEN  BY  M  OLESKIW  0N:790526  LAST  MOD  I F I  ED : 80 1 227 

C 

C  CALCULATE  TRAJECTORIES  OF  DROPLETS  IN  POTENTIAL  FLOW 
C  ABOUT  AN  AEROFOIL 
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C 

DOUBLE  PRECISION  DFLOAT . U1NF , C , RD , CD . GS . RDS , RHOA , RHOD . NUS . 
.MU,DTS(6) , DEL , XP( 7  ) ,YP(7 ) , WDSREL , DBLE , HF , UST , VST , EPS , 

. CC 1 ,CC2 ,C3 ,C4 . C5.C6.C7 . C8 . C9 . C 10 , C 1 1 ,C12 .C13.C14 .C15.C 16. 

. C 17 .C 18 .C19.C20.C2 1 ,  C22 , C23 . C24 , F , T SLOPE , HFP , ADD, 

. XL , YL . COORD .CLAP . XCDUPR . YCDUPR . YCAUPR . XCDLPR . YCDLPR . 

. YCALPR . PIM1 . PI M2 . FPIM1 , FP I M2, XCOL L . YCOLL . DABS .DSIGN, 

.LT(51 ,2) . CLAPP, K,K1 , LG , LP 1 , LTH . TOL , XN . YN , YOG . X 1 
DOUBLE  PRECISION  PSI ( 7  ) . DUADX . DVADY . DMIN1 ,DTSS , L( 5 1 ) , Y0( 5 1 ) , 

. YCDU.YCDL, XCDU , XCDL , ZZ , YCAU , YCAL , UAS( 6 ) . VAS{6) ,RED(6) , 
LEN.PRDSTI , PRDSTO , PLDST I , D I  ST , TS( 500) . YUS  1 , YLS 1 , YUS2 . YLS2 , UVAT . 

. DSORT . P INF , T INF , CR I T , XO , U SLOPE , L SLOPE , YOT( 25 , 2 ) , DDD 
DOUBLE  PRECISION  X05(6) ,UDS(6) , AN( 2,6), YDS( 6 ) , TTLACN , VPSO . NA , 

. VDS!  6 ) , HT (2.6), AO « A  1 1 A2 . BO , B 1 , B2 , B3 , E5B , CO. C 1 ,C2, 

. DM1 . DO. D 1 ,02 . E5 ,UPI . UCI , VP  I . VC I . XPI . XCI , YPI , YCI . ER 1 . ER2 , 

PRO, PLD.BETAO, YCG , TH I CK , CL APPP , SLP , YOTUX , YOTLX .LINT 
C 

REAL  XMIN, XMAX . YMIN . YMAX . SNGL . X . Y , XDSP( 1 50 ) . YDSP ( 1 50 ) . YPRE V , 

.  XPREV 

C 

INTEGER  I , CDS , XZ , YZ , IJ. IK . TRJEND , SMASH . ALMOST , AT , BOTH , ACN , 

GRAZE , IC. I  CL . ICT, ICU, IG. I U. NT. PLOT I . UX . LX , 1 1 . ICLL , 1 1 1 . 

. TRJPRA , TRJPLA , PRINT  I , PR INTO , NTRAJU , NTRAJL.TYPE.TYPE2, 

. IM4, I M3. I M2, IM1 . 10. I P 1 . IV EMP , EON . PC . INT.EO 

C 

COMMON  /EQNMN/GS , RHOA . RHOD , RDS . NUS , HF 
. /AIR/ XP, YP.DEL.PS I . TYPE2/RFL/UAS . VAS . RED . CD 
. /GRID/ XMIN, XMAX . YMIN , YMAX . XZ , YZ 
. /PV/XDS . YDS , UDS , VDS/ INTEG/AN , HT 

. /PCM/ AO. A  1 , A2.B0.B1 ,B2,83,CO,C1 .C2.DM1 ,D0.D1 ,D2, 

UP  I .UCI . VP  I . VC I . ER  1  , ER2 . XPI . XCI , YPI . YCI .UST. VST 
. /LOC/TS.DTS. I , IM4 . IM3, I M2 , I M 1 , I 0 . I P 1 

COMMON  /RKFM/CC1 , CC2 . C3 . C4 , C5 . C6 , C7 . C8 . C9 . C 10 . C 1 1.C12.C13.C14, 

. C15.C 16,C17,C18,C19,C20,C21 .C22.C23.C24 
./COL/L.YO. ICT, ICU, I CL/NOSE/ XN . YN/ SRCH/DDD , III 
C 

C  IN  TYPE  =AEROFOI L  TYPE. 

C  IN  TR JPLA=PLOT  DROPLET  TRAJECTORIES.  (0  OR  1) 

C  IN  THICK=AEROFOIL  THICKNESS  IN  %. 

C  IN  AT  =  AUTO  TRAJECTORY  MODE  (O  OR  1) 

C  IN  80TH=CALCULATE  TRAJECTORIES  TO  COLLIDE  ON  BOTH  SFCS .  (O  OR  1) 

C 

10  FORMAT (/5F6.0. 2D  10.2) 

20  FORMAT!/ 14, 2 17. 16,17, F5.0.F6.0) 

25  FORMAT! /2I 7, 13, 15, 14. 13, F6. 0,08.0. 14) 

30  FORMAT ! / 30F 10.0) 

40  FORMAT!  'OSTEP'  ,T7,  'TIME'  ,T  15,  'DTS'  ,T22,  ' XDS ' ,T31,  'YDS'  ,T40.  '  PS  I 

.  T 49 ,  ' UAS  , T58 ,  'UDS'  , T67 ,  'VAS'  ,T76.  'VDS'  . T86 .  'RED'  , T94 , 

.  ' ACCN/MOD  HIST/RHS'  ,T1  14,  'USTAB' ,T123,  ' VST AB '  ) 

50  FORMAT!'  '  . I  4 , F6 . 2 , F7 . 4 , 7F9 . 5 . F 10 . 5 , 4F9 . 5 ) 

60  FORMAT! 'OCLOSEST  APPROACH  IS  Y*'.F105.'  NO.  OF  STEPS  REQUIRED 
. 13, '  PSI*' .F8.3) 

70  FORMAT! ' 1TRAJECT0RY  STARTING  POSITION  IS  X= ' . 

F6.2,  '  Y0= '  . F  9 . 5  ) 

75  FORMAT! ' -TRAJECTORY  STARTING  POSITION  IS  X='. 

. F  6 . 2 .  *  Y0= '  . F  9 . 5  ) 

80  FORMAT!  OCOLLISION  COORDS :  X=',F10.7,#  Y='.F10.7,'  L=',F107, 

.'  NO.  OF  STEPS  REQUIRED^ ' . 13) 

90  FORMAT! 'OFIRST  TRAJECTORY  HIT  AEROFOIL') 

95  FORMAT! 'OUNEXPECTED  AEROFOIL  MISS') 

96  FORMAT!  'OYO? '  ) 

97  FORMAT! F 10.0) 

C 

C  STATEMENT  FUNCTION  TO  CALCULATE  DISTANCE  BETWEEN 
C  AEROFOIL  SLOPE  AND  TRAJECTORY 
F  (  X  )  ~ TSLOPE  *  (  X  XD  +  YL-COORD 

c 

C  INPUT  PARAMETERS 

READ! 4 .  10)UINF ,C , PINF . TINF , RD , A  1 , A2 

RE AD ( 4, 20 )COS, TRJPRA, PR  INTI  .PLOT  I  , PR  I NTO . CR I T . BE TAO 
READ! 4 , 25)NTRAJU.NTRAJL, AT .BOTH, EQN.PC.DTSS, EPS, ACN 
RE AD( 4 , 30 )X0 

C  CHECK  FOR  AUTO-TRAJECTORY  MODE 
I F  !  AT  EQ  UGOTO  700 
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1949 

1950 

1951 

1952 
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1954 

1955 

1956 

1957 

1958 

1959 

1960 
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NT  =  51 
GOTO  710 

C  CHECK  TO  SEE  IF  COLLISION  EFFICIENCIES  ARE  TO  BE  CALCULATED 
C  FOR  BOTH  SFCS . 

700  NT  =B0TH+ 1 

C 

C  NON-DIMENSIONAL  VIEWPORT  DIAGONAL  LENGTH 

710  LEN=DSQRT ( DBLE ( ( XMAX-XMIN ) **2+(YMAX-YMIN)**2) ) 

C  PRINT  LENGTH  INTERVAL  WITHIN  VIEWPORT 
PRDST I  =LEN/DFL0AT ( PR  INTI ) 

C  PLOT  LENGTH  INTERVAL  WITHIN  VIEWPORT 
PLDST I =LEN/D FLOAT ( PLOT I ) 

C  PRINT  LENGTH  INTERVAL  TO  LEFT  OF  VIEWPORT 
PRDSTO=LEN/DFLOAT ( PR  INTO ) 

C  NON-DIMENSIONAL  ACCN.  OF  GRAVITY 
GS  =  0 . DO*C/UINF/UINF 
C  NON-DIMENSIONAL  DROPLET  RADIUS 
RDS=RD* 1 .D-6/C 
DEL=RDS 
C  AIR  DENSITY 

RHOA  =  P I NF  *  1  . D3/287 .0400/ (TINF+273.  16D0) 

C  WATER  DENSITY  REF:  LIST  -  SMT 
RHOD  =  999 . 1 5D0 

C  DYNAMIC  VISCOSITY  OF  AIR  REF:  LOZOWSKI  ET  AL .  (1979) 

MU= 1 .7180-5+5. 1D-8+TINF 
C  NON-DIMENSIONAL  KINEMATIC  VISCOSITY  OF  AIR: 

NUS  =  MU/RHOA/C/U I NF 
TOL  = 1 . D - 5  *  THICK 
T YPE2  =  TYPE 
I F ( PC . NE . 2 )GOTO  420 
C 

C  DETERMINE  PARAMETERS  FOR  RUNGE -KUTTA-F EHLBERG  METHOD. 

CC1= . 25D0 
CC2=3. DO/32. DO 
C3=9. DO/32. DO 
C4  =  1932 .DO/2  197 .DO 
C5=72.D2/2197. DO 
C6=7296. DO/2197. DO 
C7=439. DO/216 .DO 
C8=8 . DO 

C9=3680. DO/513  DO 
C 10=845 . DO/ 4 104 . DO 
Cl  1=8. DO/27. DO 
C12=2 .DO 

C13=3544. DO/2565. DO 
C 14= 1859 . DO/ 4 104 . DO 
C  1 5= 1 1  .D0/40. DO 
C16=25. DO/216. DO 
C 1 7  = 1 408 . DO/2565 . DO 
C18=2197. DO/ 4 1 04 . DO 
C  19= . 2D0 

C20=16. DO/135. DO 
C21=6656. DO/12825. DO 
C22  =  2856  1  . D0/56430 . DO 
C23=9 .D0/50. DO 
C24=2. 00/55. DO 
GOTO  400 

420  I F ( PC . NE .  1  )GOTO  400 
C 

C  DETERMINE  PARAMETERS  FOR  PREDICTOR-CORRECTOR  METHOD. 

A0= 1 . DO-A 1 -A2 

B0= ( 55 . DO+9 . DO*  A  1+8 . DO* A2 )/24 . DO 

B 1  =  ( -59 .00+19  DO* A  1+32  DO* A2)/24 .DO 

B2  =  ( 37  DO- 5  DO* A  1  +  8  DO*A2)/24 .DO 

B3=( -9  DO+A 1 )/24  DO 

E5B  =  ( 25 1 .00- 1 9. DO* A  1 -8 .DO* A2)/€. DO 

C  1=A  1 

C2  =  A2 

C0=1 .D0-C1-C2 

DM  1  =  ( 9 . DO- C 1 ) /24 . DO 

D0= ( 1 9 . D0+ 13.D0*C1+8.D0*C2)/24 . DO 

D 1 = ( -5 . D0+ 13  D0*C1+32 . 00*C2 )/24 . DO 

D2=( 1 . D0-C1+8  DO*C2)/24.DO 

E 5  =  ( -  1 9 . D0+ 1 1  DO* Cl -8  DO*C2)/6.DO 
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1982 

ERt=F.5B/(E5B-E5) 

1983 

ER2=E5/(E5B-E5) 

1984 

C 

• 

1985 

C 

FOR 

EACH  TRAJECTORY  (OR  TRAJECTORY  SET): 

1986 

ENTRY  TRAJEK 

1987 

400 

IF(AT.EQ.O) GO T 0  390 

1988 

RE  ADC  4, 30)  (YO(I  Kl  =  1  .NT) 

1989 

390 

DO  200  I J= 1 , NT 

1990 

IF ( AT . EO. 1 )GOTO  395 

1991 

WR I TE ( 6 , 96  ) 

1992 

READ ( 5 , 97 ) Y0( I J) 

1993 

t F ( DABS ( Y0( IJ) )  LT .  1  .D-10)GOTO  690 

1994 

395 

I  G=  1 

1995 

GRAZE  =  1 

1996 

IC=  1 

1997 

I  NT  =0 

1998 

K  1  =  1  .  DO 

1999 

K=0 . 85DO 

2000 

YDS (  1  )  “ Y0( I J) 

2001 

C 

SET 

COUNTERS 

2002 

405 

I M4  *  2 

2003 

I  M3  =  3 

2004 

I  M2  ®  4 

2005 

I M 1  =5 

2006 

10  =  6 

2007 

IP  1  =  1 

2008 

C 

2009 

C 

DROPLET  AT  INITIAL  POSITION 

2010 

XDS ( 1 ) =X0 

201  1 

WR I T E ( 6 , 75  )  XDS(  1  )  , YDS(  1 ) 

2012 

WRITE (7, 70)  XDS (  1  )  ,YOS<  1 ) 

2013 

I F ( PC . NE . 1 )GOTO  410 

2014 

C 

2015 

SET 

PREVIOUS  PREDICTOR-CORRECTOR  VALUES  TO  0. 

2016 

XP I =0 . DO 

2017 

XC I =0 . DO 

2018 

Y  P I =0 . DO 

2019 

YC I =0 . DO 

2020 

UPI =0.00 

202  1 

UC I =0 . DO 

2022 

VP  I =0 . DO 

2023 

VC  1=0  DO 

2024 

4  10 

I  F ( ACN . EO .  1  )GOTO  415 

2025 

C 

2026 

c 

SET 

DROPLET  TRAVELLING  WITH  JUST  SLIGHTLY  GREATER  VELOCITY 

2027 

c 

THAN  AIR  ( RED=0 . 00 1 ) 

2028 

CALL  A I RVEL ( XDS (  1  )  . YDS( 1 ) ,UAS(  1  )  . VAS(  1 ) ,4 ) 

2029 

c 

CALCULATE  TOTAL  AIR  VELOCITY. 

2030 

UVA  T  =  DSQRT ( UAS (  1  )*UAS(  1  )+VAS(  1  )*VAS(  1 ) ) 

2031 

c 

CALCULATE  TOTAL  STARTING  RELATIVE  VELOCITY. 

2032 

WDSRE  L  = 1  . D-3*NUS/2 . DO/RDS 

2033 

c 

CALCULATE  INITIAL  DROPLET  VELOCITY 

2034 

UDS(  1  )  =UA  S (  1)M1. DO+WDSRE  L/UVAT ) 

2035 

VDS(  1  )=VAS(  DM1  .DO+WDSREL/UVAT  ) 

2036 

GOTO  416 

2037 

c 

SET 

GRID  FOR  INITIAL  DROPLET  VELOCITY  CALCULATIONS 

2038 

c 

2039 

4  15 

XP { 6  )  -XDS (  1  ) *2 . DO* RDS 

2040 

XP ( 7  )  =XDS (  1  ) *2 . DO*RDS 

204  1 

YP ( 6  )  =  YDS (  1  MRDS 

2042 

YP ( 7  )  = YDS (  1  ) -RDS 

204  3 

CALL  AIRVELC  XDS(  1  )  . YDS(  1  )  . UAS (  D . VAS (  1) . 7  ) 

2044 

C 

CALCULATE  DUA/DX 

2045 

DUADXM  PSI (6)*PSI( 4 )-PSl ( 7 )-PSI ( 3  D/4 . DO/RDS/RDS 

2046 

C 

CALCULATE  DVA/DY 

2047 

DVADY=( PSI ( 3)+PSI(7)*PSl(6)-PSI(4))/4  DO/RDS/RDS 

2048 

C 

TOTAL  POTENTIAL  FLOW  ACCELERATIVE  TERM 

2049 

UVA  T "DSORT ( DUADX  *  DUADX+DVADY *DVAD Y ) 

2050 

C 

CALCULATE  TOTAL  STARTING  RELATIVE  VELOCITY 

205  1 

WDSRE L=  1  . D  ~  3  *  NUS/ 2 .DO/RDS 

2052 

c 

ASSURE  STARTING  RED-0.001  WEIGHTED  BY  POTENTIAL  FLOW 

2053 

c 

ACCELERATIVE  COMPONENTS. 

2054 

UDS1  1  )=UASC  1  )  -DUADX/UVAT  *WDSREL 

2055 

VOS (  1  )  = VAS (  1  )-DVADY/UVAT  * WDSRE  L 
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2056 

416 

CALL  DRAG(UDS( 1  )  , VDS( 1  )  ,UAS( 1 ) . VASl  1  )  .CDS.RED( 1  )  .CD) 

2057 

HT ( 1 . 1 ) =0 . DO 

2058 

HT ( 2 , 1 ) =0 . DO 

2059 

C  CALCULATE  STARTING  ACCELERATIONS: 

2060 

I F ( E  ON . E  0  -  0 ) GO  T  0  417 

206  1 

EQ  =  1 

2062 

GOTO  418 

2063 

4  17 

EQ  =  0 

2064 

418 

CALL  ACCN( UDS(  1  )  , VDS( 1  )  ,UAS( 1  )  . VAS( 1 ) . RED( 1 ) , CO . EO , 0 

2065 

I F ( TR JPRA . EO. 1 ) WRITE (7, 40) 

2066 

111  =  1 

2067 

1=0 

2068 

I  K  =  0 

2069 

TR JEND=0 

2070 

Al M0ST=0 

207  1 

DT  '(  1  )  =DTSS 

2072 

CLAP= 1 . D 1 

2073 

XCDL  =0 . DO 

2074 

XCDU=0 . DO 

2075 

YCAL=0 . DO 

2076 

Y  CAU  =  0 . DO 

2077 

YCDL  =0 . DO 

2078 

Y  CDU  =  0 . DO 

2079 

SMASH=0 

2080 

T  S (  1  )=0.D0 

208  1 

PLD  =  0 . DO 

2082 

C 

2083 

100 

PRD=0 . DO 

2084 

I F ( PLD . LT . PLDSTI  )GOTO  105 

2085 

102 

PLD=0 . DO 

2086 

C 

2087 

C  INCREMENT  INDICES 

2088 

105 

I TEMP= I M4 

2089 

I M4  =  I  M3 

2090 

I  M3  =  I  M2 

2091 

I  M2  =  I M  1 

2092 

IM 1 = 10 

2093 

I0=IP1 

2094 

rP1=ITEMP 

2095 

1  =  1  +  1 

2096 

HFP=HF 

2097 

C 

2098 

C  INTEGRATE  EONS.  OF  MOTION 

2099 

I F ( PC . EO . 2 )CALL  RKF4 ( EON , CDS . EPS ) 

2  100 

IF( I  .GE .4. AND. PC. EO.  1  )CALL  PC4(EQN,CDS) 

2  101 

I F ( I . LT . 4 . AND . PC . EO . 1 . OR . PC . EO . 0 ) CA LL  RK4( EON, CDS ) 

2  102 

C 

2  103 

C  CALCULATE  DISTANCE  SINCE  LAST  PRINT/PLOT  OF  DROPLET  POSITH 

2  104 

D I  ST  =DSQRT ( ( XDS( I P 1 )-XDS( 10) )**2+( YDS( I P 1 )-YDS( 10) )* 

2  105 

PRD  =  PRD  +  D I  ST 

2  106 

X  =  SNGL ( XDS ( IP1  ) ) 

2  107 

XPREV=SNGL(XDS( 10) ) 

2  108 

I  F ( X  .GT . XM I N ) GOTO  190 

2  109 

I F ( PRD . GE . PROS TO )GOTO  230 

2  1  10 

GOTO  105 

2  111 

190 

Y  =SNGL( YDS( IP1 ) ) 

2  112 

YPREV=SNGL( YDS( 10) ) 

2  113 

C  CHECK 

FOR  OUT-OF-BOUNDS. 

2  114 

IF(Y.GE. YMAX )GOTO  2 1 1 

2  1  15 

IF(BOTH.EO. 0 ) GOTO  191 

2  1  16 

IF( Y . LT . YMIN. AND  YPREV.GT . YMIN )GOTO  212 

2  117 

19  1 

I F ( X  GE . XMAX )GOTO  213 

2  1  18 

PLD=PLD+D 1ST 

2  1  19 

IF( X  GE . SNGL ( XN )  AND . X  LE . 1  0)G0T0  240 

2  120 

I F ( IK.EO.O. AND.TRJPLA  EO . 1 )GOTO  226 

2  12  1 

I F ( PLD . GE  . PLDSTI  )GOTO  220 

2  122 

I F ( PRD . GE .PRDSTI )GOTO  230 

2  123 

GOTO  105 

2  124 

C 

2  125 

C  HOW  CLOSE  IS  DROPLET  TO  AEROFOIL? 

2  126 

C  COUNT 

NUMBER  OF  STEPS  PAST  NOSE 

2  127 

240 

ALMOST  ^ALMOST  + 1 

2  128 

IF(YDS(IP1)  LT.YN) GO T 0  310 

2  129 

IF( YDS( 10)  GT  YNJGOTO  320 
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2  130 

TSLOPE  =  ( YDS( IP1 ) -  YDS ( IO) )/(XDS( IP1 )-XDS< 10) ) 

2131 

I F ( DABS( T5L0PE -0 . DO ) . LT . 1 .D-70)GOTO  320 

2132 

X I = ( YN-YDS< 10) )/TSL0PE+XDS( IO)+RDS 

2133 

I F ( XI  GT . XN . AND . XI . LT . XDS( IP1 )+RDS/ 1 . D6 )GOTO  330 

2  134 

C 

2135 

C  FOR  UPPER  SFC . : 

2  136 

320 

CALL  SFC( XDS( IP  1  )  , YUS  1 , 1 .0, ZZ ) 

2  137 

CALL  SFC ( XDS ( IP1 ) +RDS . YUS2 . 1 .0,22) 

2  138 

USLOPE=DSQRT (RDS*RDS+( YUS2-YUS  1  )*  +  2) 

2  139 

C  PREVIOUS  CLOSEST  APPROACH  X  AND  Y  COORDS. 

2  140 

XCDUPR=XCDU 

2  14  1 

YCDUPR-YCDU 

2  142 

YCAUPR=YCAU 

2  143 

C  CALCULATE  DROPLET  Y  COORD.  OF  CLOSEST  APPROACH 

2  144 

Y  CDU  =  YDS ( IP1 ) -RDS *RDS/USLOPE 

2  145 

C  CALCULATE  DROPLET  X  COORD.  OF  CLOSEST  APPROACH 

2  146 

XCDU  =  XDS ( I P 1  )  +  RDS  * ( YUS2-YUS1 )/USLOPE 

2  147 

C  CALCULATE  AEROFOIL  X  AND  Y  COORDS.  OF  CLOSEST  APPROACH 

2148 

CALL  SFC (XCDU.YCAU, 1,0, ZZ) 

2  149 

C  STORE 

CLOSEST  APPROACH  VALUE 

2  150 

CL AP  =  DM I N 1 ( CLAP . ( YCDU-YCAU)  ) 

2  15  1 

C  CHECK 

FOR  DROPLET-AEROFOIL  'COLLISION' 

2  152 

I Fi ( YCDU-YCAU ) . LE . RDS*CRIT/ 1 .D2 )GOTO  505 

2  153 

I F ( PLD  GE . PLDSTI  )GOTO  220 

2  154 

I F ( PRD .  GE  . PRDST I  )GOTO  230 

2  155 

GOTO  105 

2  156 

C 

2  157 

C  COLLISION  FLAGGED 

2  158 

505 

$MASH=  1 

2  159 

IF( YCDU.GT. YCAU)GOTO  520 

2  160 

IF( ALMOST  EO. 1 )GOTO  500 

2  16  1 

XL  =  XCDUPR 

2  162 

Y  L  =  YCDUPR 

2  163 

TSLOPEM  YCDU-YCDUPR )/ ( XCDU- XCDUPR ) 

2  164 

P I  M2  - XDS ( 10) 

2  165 

CALL  SFC(PIM2. COORD, 1 ,O.ZZ) 

2  166 

F  P I  M2  =  F ( P I  M2  ) 

2  167 

P I M 1 =XDS ( I 0 ) +RDS 

2  168 

CALL  SFCIPIM1 , COORD, 1 ,0,ZZ) 

2  169 

FPIM1 =F ( PIM 1 ) 

2  170 

GOTO  510 

2  17  1 

C  NEAR 

NOSE  COLLISION 

2  172 

500 

XL -XDS ( 10 ) +RDS 

2  173 

YL  =  YDS (  10) 

2  174 

T  SLOPE  =  ( YDS (  IP1  )  -  Y  L  )  / ( XDS (  IP1  )+RDS-XL) 

2  175 

IF (DABS( TSLOPE-O. DO )  LT . 1 .D-7O)G0T0  507 

2  176 

P  I  M2  =  XN 

2  177 

COORD  =  YN 

2  178 

FPIM2=F(PIM2 ) 

2  179 

PIM1-XDS(  IP  1  ) 

2  180 

COORD  =  YUS  1 

2  18  1 

FPIM1=F(PIM1  ) 

2  182 

GOTO  510 

2183 

507 

XCOLL  =  XN 

2  184 

YCOLL  =  YN 

2  185 

LTH=0. DO 

2  186 

GOTO  210 

2  187 

C  AN  'ALMOST'  COLLISION 

2  188 

520 

XCOLL  =  XDS(  IP  1  ) 

2  189 

CALL  SFC( XCOLL , YCOLL .  1  ,  1 , LTH) 

2  190 

GOTO  210 

2  19  1 

C 

2192 

C  ITERATE  TO  COLLISION  LOCATION  USING  SECANT  METHOD. 

2  193 

5  10 

XCOLL  r  P I M 1 -FPIM1 ♦(PIM1-PIM2 )/( FPIM1 -FPIM2) 

2  194 

I F ( XCOLL . GT . XN )GOTO  511 

2  195 

XCOLL  =  XN 

2  196 

COORD= YN 

2  197 

GOTO  512 

2  198 

5  1  1 

CALL  SFC(XCOLL. COORD. 1 .O.ZZ) 

2  199 

512 

PIM2-PIM1 

2200 

F  P  I  M2  =  F  P  I  M  1 

2201 

PIM1-XC0LL 

2202 

FPIM1~F(XC0Lt  ) 

2203 

I F ( DABS ( F  P I M 1  )  GT  RDS * CR I T /  1  D2 ) GO TO  510 

2204 

CALL  SFC( XCOLL , YCOLL , 1 . 1 , LTH ) 

2205 

GOTO  210 

2206 

310 

IF(YOS(IO)  .LT.YN1G0T0  330 

2207 

TSLOPE  =  ( YD$( IP1  )  -  YDS ( 10) )/(XDS( IP1 )-XDS( 10) ) 

2208 

IF(DABS(TSL0PE-0. DO ) . LT . 1 . D- 70) GOTO  330 

2209 

XI  =  ( YN-YDS( IO) ) /TSLOPE  +  XDS ( IO)+RD$ 

22  10 

I F  < XI  . GT . XN . AND . XI . LT . XDS( IP1  )  +  RDS/  1  . D6)G0T0  320 

22  1  1 

C 

22  IP- 

C  FOR 

LOWER  SFC . : 

2213 

330 

CALL  SFC(XDS( IP1 ) , YLS1 ,0.0, ZZ ) 

2214 

CALL  SFC ( XDS ( I P 1 )+RDS,YLS2,0,0. ZZ) 

2215 

LSLOPE=DSQRT( RDS*RDS+( YLS2-YLS  1  )**2  ) 

22  16 

C  PREV 

IOUS  CLOSEST  APPROACH  X  ANO  Y  COORDS 

22  17 

XCDLPR=XCDL 

22  18 

YCDLPR=YCDL 

22  19 

YCALPR- YCAL 

2220 

C  CALCULATE  DROPLET  Y  COORD.  OF  CLOSEST  APPROACH 

222  1 

YCDL  =  YDS ( IP  1  )  +  RDS*RDS/LSLOPE 

2222 

C  CALCULATE  DROPLET  X  COORD.  OF  CLOSEST  APPROACH 

2223 

XCDL  =  XDS ( IP1 )  +  RDS*(YLS1-YLS2 )/LSLOPE 

2224 

C  CALCULATE  AEROFOIL  X  AND  Y  COORDS.  OF  CLOSEST  APPROACH 

2225 

CALL  SFC ( XCDL .YCAL.O.O.ZZ) 

2226 

C  STORE  CLOSEST  APPROACH  VALUE 

2227 

CLAP=DS IGN( CLAP , -1.00) 

2228 

CLAP -DMAX 1 (CLAP , ( YCDL -YCAL  )  ) 

2229 

C  CHECK  FOR  DROPLET-AEROFOIL  ' COLLISION ' 

2230 

IF( (YCAL -YCDL ) . LE . RDS  *CR I T/ 1  .D2)G0T0  504 

223  1 

IF(PLD.GE.PLDSTI  )GOTO  220 

2232 

I F( PRD . GE . PRDSTI )GOTO  230 

2233 

GOTO  105 

2234 

C 

2235 

C  COLLISION  FLAGGED 

2236 

504 

SMASH= 1 

2237 

IF( YCAL . GT . YCDL )GOTO  570 

2238 

IF ( ALMOST . EO . 1)G0T0  550 

2239 

XL  =XCDLPR 

2240 

YL=YCDLPR 

224  1 

TSLOPE  tfCDL- Y  CDLPR ) / { XCDL -XCDL PR ) 

2242 

P I M2  =  XDS ( 10) 

2243 

CALL  SFC( PIM2. COORD, 0,0, ZZ ) 

2244 

F  P I  M2  =  F ( P I  M2 ) 

2245 

P I M 1 =XDS ( 10 ) +RDS 

2246 

CALL  SFC(PIM1 .COORD, 0.0, ZZ) 

2247 

FPIM 1 =F( PIM1 ) 

2248 

GOTO  560 

2249 

C  NEAR 

NOSE  COLLISION 

2250 

550 

XL=XDS( 10 ) +RDS 

225  1 

YL  =  YDS(  10) 

2252 

TSLOPE  =  ( YDS ( IP1 ) - YL ) / ( XDS ( I P 1 )+R0S*XL) 

2253 

I F ( DABS ( TSLOPE -0 . DO) . LT . 1 . D -70) GOTO  556 

2254 

P I  M2 *  XN 

2255 

COORD* YN 

2256 

F P I M2 “F ( P I M2  ) 

2257 

P IM 1 =XDS ( I P  1  ) 

2258 

COORD  =  YLS 1 

2259 

FPIM 1 =F ( PIM1 ) 

2260 

GOTO  560 

226  1 

556 

XCOL  L  =  XN 

2262 

YCOLL  =  YN 

2263 

L  TH  =  0 . DO 

2264 

GOTO  210 

2265 

C  AN  ' ALMOST '  COLLISION 

2266 

570 

XCOLL “ XDS ( IP1  ) 

2267 

CALL  SFC(XCOLL. YCOLL. 0. 1 .LTH) 

2268 

GOTO  210 

2269 

c 

2270 

C  ITERATE  TO  COLLISION  LOCATION  USING  SECANT  METHOD 

227  1 

560 

XC0LL  =  PIM1 -FPIM1 M  PI  Ml -PI  M2 )/( FPIM1 -FPIM2 ) 

2272 

IF( XCOLL .GT . XN  )GOTQ  561 

2273 

XCOLL - XN 

2274 

COORD  =  YN 

2275 

GOTO  562 

2276 

56  1 

CALL  SFC( XCOLL .COORD .0.0. ZZ ) 

2277 

562 

P I M2  =  P I M 1 

61 


2278 

2279 

2280 
228  1 
2282 

2283 

2284 

2285 

2286 

2287 

2288 

2289 

2290 

2291 

2292 

2293 

2294 

2295 

2296 

2297 

2298 

2299 

2300 

2301 

2302 

2303 

2304 
2  305 
2  306 

2307 

2308 
2  309 
2310 
23  1  1 
23  12 

2313 

2314 

2315 

2316 
23  17 
23  18 
23  19 
2320 

232  1 

2322 

2323 

2324 

2325 

2326 

2327 

2328 

2329 

2330 

233  1 

2332 

2333 

2334 

2335 

2336 
2  3  37 

2338 

2339 

2340 

234  1 

2342 

2343 

2344 

2345 

2346 

2347 

2348 

2349 

2350 
2  35  1 


FPIM2=FPIM1 
P I M 1 -XC0LL 
FPIM1=F( XCOLL ) 

I F ( DABS ( FP I M 1  )  . GT . RDS  *CR I T/ 1 .D2)G0T0  560 
CALL  SFC(XCOLL, YCOLL.O, 1 . LTH) 

C 

C  END  OF  TRAJECTORY  FLAGGED:  COLLISION 

210  TR JEND= 1 

I F( IK . EQ.O.OR . TRJPLA . EO  Q)GOTO  230 
IK  =  IK+  1 

XDSP ( IK ) -$NGL( XCOLL ) 

YDSP ( IK ) =SNGL ( YCOLL ) 

GOTO  230 

C  END  OF  TRAJECTORY  FLAGGED:  EXCEEDED  YMAX 

211  TR JEND= 1 

IF( IK . EO .0 .OR . TRJPLA . EQ.O)GOTO  230 
I K  = I K  + 1 

XDSP ( IK)=(X-XPREV)/( Y-YPREV)*( YMAX - YPRE V ) +XPRE V 
YDSP ( IK  )  = YMAX 
GOTO  230 

C  END  OF  TRAJECTORY  FLAGGED:  EXCEEDED  YMIN 

212  TR JEND= 1 

I F ( IK. EQ.O.OR. TRJPLA . EO-0)GOTO  230 
IK=IK+1 

XDSP( IK)=(X-XPREV)/( Y-YPREV)*( YMIN- YPREV ) +XPREV 
YDSP ( IK)  =  YM I N 
GOTO  230 

C  END  OF  TRAJECTORY  FLAGGED;  EXCEEDED  XMAX 

213  TR JEND  = 1 

IF ( IK. EQ.O.OR . TRJPLA . EQ.O)GOTO  230 

I K  = I K  + 1 

XDSP ( IK) =  XMAX 

YDSP ( IK )  =  ( Y-YPREV ) / ( X-XPRE V ) * ( XMAX -XPRE V ) +YPRE V 
GOTO  230 
C 

C  STORE  PLOT  COORDINATES  FOR  FIRST  POINT  WITHIN  WINDOW 
226  I K  = I K+  1 

XDSP ( IK ) =XMIN 

YDSP( IK  )  =  ( Y-YPREV)/ ( X -XPREV ) * ( XMIN-XPREV )+YPREV 
GOTO  230 

C  STORE  COORDS  FOR  LATER  PLOTTING 
220  IF(TRJPLA.EO.O)GOTO  230 

I K  =  I K  +  1 

XDSP ( IK) =SNGL ( XDS ( 10) ) 

YDSP ( I K ) =  SNGL ( YDS ( 10)) 

230  I F ( TR JPRA . EO . 0 . AND . TR JEND . EO . 0)G0TD  100 

IF ( PRD . LT  .  PRDSTI . AND . TR JEND . EQ.O)GOTO  102 

C 

C  PRINT  INTERVAL  EXCEEDED 

TTL ACN=DSQRT ( AN( 1 , IO)*AN{  1 . 10  )  +  AN( 2 . 1 0 )  +  AN( 2 . 10 ) ) 
VPSQ=UOS( 10 ) *  UDS  ( I O )  +  VDS( IO)+VDS( 10) 

NA=RDS*TTLACN/DTS( IO)/VPSO 
I F ( TRJPRA . EO.O)GOTO  181 
C 

C  PRINT  TRAJECTORY  INFO 

I F ( PC . EG. 1 . AND . I .GT . 4 )GOTO  235 

WR I T  E ( 7 , 50 )I , TS( I ) ,DTS( IO),XDS( 10) ,YDS( 10) ,PSI(5) ,UAS( 10) 
UDS ( 10) , VAS( 10) , VDS ( 10) . RED( IO) , NA . HFP 
I  F  { TRJEND . EQO )GOTO  100 
GOTO  225 

235  WRITE (7. 50) I ,  TS( I ).DTS( 10 ) , XDS ( 10 ) , YDS ( 1 0 ) , PS  I ( 5 ) , UAS ( 10) 

UDS ( 10) . VAS( 10) ,VDS( 10) ,REO( I O ) . NA . HFP , UST . VST 
I F( TRJEND  EOO)GOTO  100 
225  1=1+1 

WRI TF( 7 . 50) I  . TS( I  )  ,DTS(  IP  1  )  . X . Y 
181  I F ( TRJPLA . EQ.O)GOTO  180 

C 

C  PLOT  TRAJECTORIES: 

XDSP (  IK+ 1 )  =  XM I N 

XDSP ( IK+2 )=( XMAX-XMIN)/2 1 .0 

YDSP ( I K  +  1  )  =  YMI N 

YDSP (  IK  +  2  )  =  ( YMAX -YMIN)/ 10 . 5 

CALL  LINE1XDSP. YDSP. IK. 1.0,0) 

180  I F ( SMASH . EG . 1 )GOTO  195 


62 


2352 

2353 

2354 

2355 

2356 

2357 

2358 

2359 

2360 

2361 

2362 

2363 

2364 

2365 

2366 

2367 

2368 

2369 

2370 

237  1 

2372 

2373 

2374 

2375 

2376 

2377 

2378 

2379 

2380 

238  1 

2382 

2383 

2384 

2385 

2386 

2387 

2388 

2389 

2390 

2391 

2392 

2393 

2394 

2395 

2396 
2  397 

2398 

2399 

2400 

2401 
2  402 

2403 

2404 
2  405 
24  36 
2407 
2403 
2409 
24  10 
24  1  1 
24  12 
24  1  3 
24  14 
24  15 
24  16 
24  17 
24  18 
24  19 
24  20 
242  1 

2422 

2423 

2424 
2  4  2  ^ 


WR I T  E ( 6 , 60 )CL AP , I . PSI ( 5  ) 

WR I T  E ( 7 , 60 )CL AP , I , PSI (5) 

GOTO  196 

195  WRITE(6 , 80 )X COLL , YCOLL , LTH, I 
WRITE( 7 , 80 )X COLL , YCOLL , LTH. I 

196  IF(AT.EOO) GOTO  200 

I F ( GRAZE . EQ - O )GOTO  630 
IF( SMASH. EQ, 1 )GOTO  610 
C 

C  ITERATE  TOWARD  THE  GRAZING  TRAJECTORY 
IF(  IG. EO.  1  ) GO  T  0  600 
IF (DABS (CLAP)  . LE . TOL  )K  =  K 1 

C  FIND  FIND  NEW  YO  POSITION  BY  USING  THE  SECANT  METHOD  TO  ESTIMATE  THE 
C  LOCATION  OF  YO  AT  GRAZING 

SLP  =  ( YOT ( IG. I J ) - YOT ( IG~ 1 , I J ) )/ (CLAP -CLAPP ) 

I F ( DABS ( SLP  )  .  LT  1 .200. OR.  IG.LT.3)G0T0  340 
SLP  =  ( YOT ( IG. I J ) -YOT ( IG-2 . Id ) )/ (CLAP-CLAPPP ) 

K  =  K  1 

340  YOT (  IG+ 1 , I J  )  =  YOT ( I G  t I J ) -K+CLAP *SLP 

C  SET  PREVIOUS  CLOSEST  APPROACH 
CL APPP=CLAPP 
CLAPP  =  CLAP 
IG= IG+ 1 

YDS ( 1 )=YOT( IG, I J) 

GOTO  405 

C  AFTER  FIRST  MISSING  TRAJECTORY,  ESTIMATE  NEW  YO  VIA  CLAP 
600  YOT (  1 . I J  )  =  YO ( IJ) 

I F ( DABS ( CLAP ) . LE .TOL )K=K 1 
YOT ( 2 , I J )  =  YOT (  1  . I J  )  -K*CL AP 
C  L APP  =  CL AP 
I  G  =  2 

YDS ( 1 )=Y0T(2, IJ) 

GOTO  405 

C 

C  THIS  IS  THE  GRAZING  TRAJECTORY 
610  I F ( IG.GT  1 ) GOTO  625 

WRITE  (O’.  90) 

C  ADJUST  F'RST  TRAJECTORY  TO  BE  A  NEAR  MISS 
I F ( I J . EO  2 ) GOTO  605 
YO (  1 )  =  Y0(  1  )  +5 . D  -  4 
GOTO  606 

605  Y0(  1  )  = Y0(  1  )  -5 . D 7 4 

606  YOS (  1  )  =  YO (  1  ) 

GOTO  405 

625  GRA  ZE  =0 

YOG  =  YOT ( IG .  I J ) 

YOT (  1 ,IJ)  =  YOG 
YCG= YCOLL 
LG=LTH 
LP  1  =  LG 


C 

C  THESE 
630 


635 


8  10 


ARE  COLLIDING  TRAJECTORIES 
IF( SMASH. EQ.  1  )GOTO  635 
WR I T  E ( 6 , 95 ) 

IC= IC- 1 
GOTO  640 

IF( IC . EO.  1  ) GOTO  800 

I F ( DABS ( DS I GN( YCOLL ~YN . YCG- YN ) -YCOLL +YN )  GT . 1 . D- 10) GOTO  645 
I F (  I  NT . EQ . 0  )GOTO  8  10 
I  NT  =0 

L  T (  IC.  I J  )  =  L  TH 
IC  =  IC+  1 
GOTO  820 

IF(LT(IC-1.IJ)-LTH.LE.1. 35D0*LINT )GOTO  800 

LT(  IC+  1  ,  I J  )  =  L  T  H 

YOT  (  IC+1  I  J  )  =  YOT  (  IC,  IJ) 

I  NT  -  1 

YOT (  IC .  I J ) "0 . 6D0*  YOT ( IC“ 1 .  I J ) +0 . 400*  YOT (  IC .  I J) 

YDS ( 1 )=Y0T( IC. IJ) 

GOTO  405 
L  T ( IC.  I J )  =  L  T  H 
IM  IC  GT  1  ) GO T 0  633 
I F ( BOTH  EQ . O ) GOT 0  63  1 
IF ( I J . EQ . 2 )GOTO  632 


800 

820 


2426 

C  ESTIMATED  INTERVAL  IN  L  BETWEEN  COLLISIONS. 

2427 

LINT=LG/(DFLOAT( NTRA JU ) + 1 .DO) 

2428 

ADD* -0 . 5D0 

2429 

GOTO  633 

2430 

632 

LINT=LG/(DFL0AT(NTRAJL)+1 .DO) 

243  1 

ADD  * -0 . 5D0 

2432 

GOTO  633 

2433 

63  1 

LINT=LG/(DFLOAT(NTRAJU )+0. 5D0) 

2434 

ADD  =  0 . 5D0 

2435 

633 

LP 1  *  LP 1  - L I NT 

2436 

IF ( IC . EO . 1 )LP1=LP 1 -ADD* LI NT 

2437 

I F ( DABS ( DS I GN( LP 1 , LTH ) -LP 1 ) . GT . 1 . D- 10. OR . DABS( LP 1 ) . LT 

2438 

GOTO  640 

2439 

IC=IC+  1 

2440 

C 

244  1 

C  ESTIMATE  NEW  VO  TO  SPREAD  POINTS  EVENLY  ALONG  CE  VS  L  CURVE 

2442 

IF  (BOTH. EO. 1 )GOTO  620 

2443 

YOT( IC. I J)  =  2 . DO*  YOG/ LG* LP 1  * (  1  . DO- LP 1 / 2 . DO/LG ) 

2444 

YDS ( 1 ) = Y01 ( IC. I J) 

2445 

GOTO  405 

2446 

620 

IF( IJ.EQ.2  )GOTO  850 

2447 

YOT ( IC.  1 ) =  BE  T  AO*  LP 1  * (  1  . DO- LP 1 /2  DO/ LG ) 

2448 

+ YOG -BET AO* LG/ 2 . DO 

2449 

GOTO  860 

2450 

850 

YOT (  I C . 2 ) * -BET  AO*  LP 1  * (  1  . DO-LP 1 /2 . DO/LG ) 

245  1 

+  YOG+BE  T  AO*  LG/ 2 . DO 

2452 

IF ( YOT ( IC . 2  )  LT . YOT( ICU. 1 ) )GOTO  860 

2453 

IC=IC- 1 

2454 

GOTO  640 

2455 

860 

YDS (  1  )=Y0T( IC. IJ) 

2456 

GOTO  405 

2457 

640 

IF( IJ.EO. 1 )ICU=IC 

2458 

I  F ( IJ.E0.2)ICL  =  IC 

2459 

GOTO  200 

2460 

645 

L  T (  IC, I J  )  =  LTH 

246  1 

I F ( I J . EO . 2 ) GOTO  646 

2462 

ICU=IC- t 

2463 

UX=  1 

2464 

GOTO  200 

2465 

646 

I  CL  =  I C -  1 

2466 

LX  =  1 

2467 

200 

CONTINUE 

2468 

C 

2469 

C  TRANSFER  COLLISION  INFO  TO  SINGLE  MONOTONI C AL L Y  INCREASING 

2470 

C 

( IN  L )  VECTORS 

247  1 

I F ( BOTH . EO . 1 ) GOTO  660 

2472 

I F ( DABS ( LT ( ICU , 1 ) -0 .DO ) . GT . 1  D-4)G0T0  651 

2473 

ICU=ICU- 1 

2474 

651 

Y0( ICU  + 1 ) =0 . DO 

2475 

L  (  ICU+1  )  =0  .  DO 

2476 

DO  650  1  =  1  . ICU 

2477 

IU  =  2  * ICU+2- I 

2478 

Y0( I ) = - YOT ( 1.1) 

2479 

Y  0 (  I U )  =  YOT (  1.1) 

2480 

L ( I  >  =  -LT ( I ,  1  ) 

248  1 

L ( I U  )  =  LT (  1.1) 

2482 

650 

CONTINUE 

2483 

I CT  * 2  *  I CU  + 1 

2484 

I  CL  =  ICU*- 1 

2485 

GOTO  690 

2486 

660 

IF ( UX  EO .  1 ) YOTUX  *  YOT ( ICU  + 1.1) 

2487 

I F ( LX . EO .  1  )  YOTLX  =  YOT( ICL  +  1  , 2  ) 

2488 

11=0 

2489 

DO  670  1  =  1.  I CL 

2490 

I F ( UX  NE . 1 )GOTO  665 

249  1 

IF< YOTUX . GE . Y0T( I  , 2 )  )GOTO  665 

2492 

IF (DABS ( YOTUX- YOT ( I . 2 ) > . LT . 1 . D-5 )GOTO  666 

2493 

IF ( DABS (YOTUX -YOT ( 1-1. 2)). LT.  t.D-5  )GOTO  666 

2494 

11=11+1 

2495 

Y0(  I I  )  =  YOTUX 

2496 

L( I I ) = - LT ( ICU  + 1.1) 

2497 

666 

UX=0 

2498 

665 

11=11+1 

2499 

Y0( I I  )  =  YOT (  1.2) 

D-  10) 


64 


2500 

2501 

2502 

2503 

2504 

2505 

2506 

2507 

2508 
2503 
2510 
25  1  1 

2512 

2513 
25  14 

2515 

2516 

2517 
25  18 
25  19 
2520 
252  1 

2522 

2523 

2524 

2525 

2526 

2527 

2528 

2529 

2530 

2531 

2532 

2533 

2534 

2535 

2536 

2537 

2538 

2539 

2540 

254  1 

2542 

2543 

2544 

2545 

2546 

2547 

2548 

2549 

2550 

255  1 

2552 

2553 

2554 

2555 

2556 

2557 

2558 

2559 

2560 

256  1 

2562 

2563 

2564 

2565 

2566 

2567 

2568 

2569 

2570 

2571 

2572 

2573 


L(II)--LT(I,2) 

670  CONTINUE 

I F ( UX . NE .  1  ) GOTO  667 

IF(DABS( YOTUX-YOT ( ICL  f  2 ) ) .  LT  ,  1  D-5)G0T0  667 

11=11+1 

Y0(  II  )  = YOTUX 

L ( 1 1 )  =  -LT ( ICU+ 1 , 1) 

667  ICLL=ICL 
ICL= I  I 

DO  680  IM,  ICU 

IU=ICU+1-I 

I F ( L X  NE .  1 1 GOT 0  675 

IF( YOTLX  GE  Y0T( IU. 1 ) )GOTO  675 

I F ( DABS ( YOTLX - YOT ( IU,  1  )  )  . LT .  1 . D-5)G0T0  676 

IF( DABS ( YOTLX -YOT( IU+ 1 . 1 ) ) . LT . 1 .D-5  )G0T0  676 

11=11+1 

Y0( II ) = YOTLX 

L (  I  I  )  =LT  ( ICLL+1 ,2) 

676  LX=0 

675  11=11+1 

Y0(  I  I  ) =  YOT ( IU.  1  ) 

L (  I  I )  =  L  T  ( IU.  1  ) 

680  CONTINUE 

I CT  =  I  I 
ICU=ICT-ICL 
690  RETURN 
END 
C 
C 

SUBROUTINE  ACCN( UD . VD . UA , VA , RED . CD , EON . T . G ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:  801216  LAST  MOD I F I  ED : 80 1 2 2 3 

C 

C  CALCULATES  RHS  OF  NON-DIMENSIONAL  EONS.  OF  MOTION 

C 

DOUBLE  PRECISION  RE L VE L , RED , NUS , RDS , APU , APV , BPU , BPV 
. ,AN( 2,6) ,  HF  , HX , HY , HT ( 2 , 6 ) . DSQRT , AU , A V . BU , BV , RHOA , 

. RHOD.GS, ALPHAR , PI , CD , UD , VD . UA , VA , TS ( 500 ) . DTS ( 6 ) . T 

C 

INTEGER  EON.G, I . IM4 , I M3 , IM2 . IM 1 , 10 . I P 1 

C 

COMMON  ALPHAR , P I /EQNMN/GS . RHOA . RHOD . RDS . NUS . HF 
. /INTEG/AN.HT/LOC/TS.DTS. I , IM4 , IM3 . IM2 . IM I . 10. IP  1 

C 

C  IN  UD  = 

C  IN  VD-DROPLET  VELOCITY  COMPONENTS. 

C  IN  UA  = 

C  IN  VA  =  A  I R  VELOCITY  COMPONENTS. 

C  IN  RED=RELATIVE  MOTION  REYNOLDS  NO. 

C  IN  CD=DRAG  COEFFICIENT. 

C  IN  EQN=PARAMETER  TO  DETERMINE  TERMS  USED  IN  EON.  OF  MOTION. 

C  IN  T  =  T I  ME  AT  THIS  TIME  STEP. 

C  IN  G=0: EXTRAPOLATE  HISTORY  TERM  SEQUENCE. 

C  IN  1 : CALCULATE  NEW  HISTORY  TERM  VALUE. 

C 

RELVEL=RED*NUS/RDS/2 . DO 
I F ( EON. E0.0)G0T0  100 

C 

C  FIRST  TWO  TERMS  IN  EON  OF  MOTION  INCLUDING  GRAVITATION  AND 
C  STEADY  STATE  DRAG.  (INCLUDES  BUOYANCY  AND  INDUCED  MASS  EFFECTS) 
APU  =  2  DO* (RHOD -RHOA )/( 2  DO*RHOD  +  RHOA ) *GS *DS I N( ALPHAR ) 

APV  =  2 . DO* ( RHOD- RHOA )/( 2 . DO*  RHOD  +  RHOA ) *GS*DCOS ( ALPHAR ) 

BPU  =  0 . 7  5DO*CD  *  RHOA /RDS/ ( 2 . DO*  RHOD  +  RHOA ) 

.  MUD-UA  )*RELVEL 

BPV  =0 . 75DO*CD*RHOA/RDS/( 2  DO*RHOD  +  RHOA ) 

.  *( VD-VA  )*RELVEL 
AN(  1  .  IP1  )=APU-BPU 
AN( 2 .  IP  1  )  =  - APV -BPV 
IF (EON  EQ  2)G0T0  300 
HF  =0 . DO 
RETURN 
C 

C  THIRD  (HISTORY)  TERM  FOR  SHEDDING  OF  VORTICITY 
300  CALL  HIST(T.G) 
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HX=-9.DO*RHOA/(2.DO*RHOD+RHOA)/RDS*DSQRT(NUS/PI )*HT( 1 . IP1 ) 

HY  =  -9 . D0*RH0A/ ( 2 . DO*RHOD+RHOA )/RDS*DSQRT (NUS/PI ) *HT ( 2 , IP  1 ) 

AN( 1 , IP1 ) =  AN( 1 . IP1 )+HX 
AN( 2 , I P 1 ) =  AN( 2 , I P 1 )+HY 
I F ( G . EQ . 0 ) RE TURN 

HF  =DSQRT ( (HX*HX+HY*HY )/( ( APU-BPU ) * *2+( APV+BPV ) ♦ *2 ) ) 

RETURN 

C 

C  FIRST  TWO  TERMS  IN  EON.  OF  MOTION  WITHOUT  BUOYANCY  AND  INDUCED  MASS 
100  AU=GS*DSIN( ALPHAR ) 

AV=GS*DCOS( ALPHAR ) 

BU=0. 375DO*RHOA/RHOD*CD/RDS*(UD-UA)*RELVEL 
BV  =  0. 375DO*RHOA/RHOD*CD/RDS*( VO-VA )*RELVEL 
AN( 1 . IP1 )=AU-BU 
AN ( 2 , IP1 ) - -AV-BV 
HF  =  0 . DO 
RETURN 
END 
C 
C 

SUBROUT  I NE  A I RVEL ( X , Y , UA5 , V AS , NP  ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  0N:8OO222  LAST  MOD I F I  ED : 80 1 2 1 6 
C 

C  CALCULATES  THE  AIR  VELOCITY  COMPONENTS  AT  A  GIVEN  LOCATION 

C 

DOUBLE  PRECISION  X , Y , UAS , VAS , XP ( 7  )  , YP ( 7  )  , XC (  10 1 ) , YC (  1 0 1 ) 

. .DEL, GAMMA ( 101 ) .D( 100),K( 101 ) ,  PI , PJKA . PUKE , 

. S I ( 100) ,C0( 100) ,PSI (7) ,DXC,DYC,DELTA,A,B,R1S,R2S, 

. R3S .DAT AN, T3, DABS, DS I GN , ALPHAR. T1 . T 2 . DLOG , R . DCOS . DS I N 

C 

INTEGER  L , NP , J , NCOU , NCOL , N , TYPE 

C 

COMMON  ALPHAR . P I / AERO3/NC0U . NCOL/ AER02/XC . YC . GAMMA , D . S I , CO 
. / A  I R/XP . YP , DEL , PSI , TYPE 

C 

C  I N  X  = 

C  IN  Y  =COORDS .  AT  WHICH  AIR  VELOCITY  15  TO  BE  DETERMINED. 

C  OUT  UAS  = 

C  OUT  VAS  =  COMPONENT  $  OF  AIR  VELOCITY. 

C  IN  NP=NUMBER  OF  POINTS  AT  WHICH  TO  CALCULATE  PSI. 

C 

N = NCOU + NCOL - 2 

C  SET  GRID  FOR  AIR  VELOCITY  CALCULATIONS 
XP(  1  )=X+DEL 
XP( 2)=X-DEL 
XP( 3  )=X 
XP( 4  )  =X 
X  P  (  5  )  =  X 
YP(  1  )  - Y 

Y  P  (  2  )  =  Y 

Y P ( 3  )  = Y  +  DE  L 
yp(4)=Y-DEL 

Y  P  (  5  )  =  Y 

DO  110  J=1 .NP 
IF ( TYPE . EQ .  -  1 )GOTO  1 15 
PSI ( J )=0.0 
DO  120  L=1.N 

C  FIND  DISTANCE  BETWEEN  CONTROL  PT .  L  AND  GRID  PT .  I.J. 

DXC=XP( J)-XC(L) 

DYC=YP(J)-YC(L) 

C  CALCULATE  COMPONENTS  OF  EQN .  9  AND  FIG.  2 
DELTA=D( L  )/2 . DO 
B  =  DXC*CO( L  )+DYC*SI ( L  ) 

A  =  D YC  * C0( L ) -DXC*SI ( L ) 

R 1S=A*A+(B+DELTA)*(B+DELTA ) 

R2S  =  A*A+(B -DELTA  )  MB -DELTA  ) 

R3S=A*A+B*B~DELTA*DELTA 
IF(R3S.LT. 1.D- 30 ) GO  TO  130 
T3  =  DA  T  AN( 2 . DO  * A *DE L TA/R3S ) 

GO  TO  140 

130  I F ( DABS( A  )  . LT .  1 . 0-30 )G0  TO  150 

T 3 -DAT AN (  ( B  +  DE L T A ) /A ) - DAT AN( ( B -DELTA )/A  ) 

GO  TO  140 
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150  T3=DSIGN( PI  ,  A  ) 

140  T  1MB+DELTA)  *DL0G(  R  IS ) 

T2=(B-DELTA) *DL0G( R2S ) 

K( L ) * ( T 1-T2+2. D0*A*T3-4 . D0*DELTA )/4 . 00/ PI 
PSI(J)"PSI(J) -GAMMA ( L ) *K( L ) 

120  CONTINUE 

R=YP(U)*DCOS( ALPHAR )-XP(J)*DSIN( ALPHAR ) 

C  ASSURE  THAT  P$I  ON  AEROFOIL  =  0. 

PS  I ( J)=PSI( J)  +  R-GAMMA(N+1 ) 

GOTO  110 

1 15  PSI ( J)=YP( J)-YP( J)/4 .D0/( (XP( J)-5 .D- 1 )**2+YP(U) *YP(U ) ) 

110  CONTINUE 

C 

C  CALCULATE  AIRSPEED  FROM  STREAMFN. 

UAS=(PSI(3)-PSI(4))/2. DO/DEL 
VAS=(PSI(2)-PSI( 1 ) )/2. DO/DEL 
RETURN 
END 
C 
C 

SUBROUTINE  DRAG( UDS . VOS , UAS . VAS . CDS , RED . CD ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:800222  LAST  MOD I F I  ED : 80 1 2 1 6 
C 

C  CALCULATES  THE  REYNOLDS  NUMBER  AND  DRAG  COEFFICIENT  OF  THE  DROPLET 

C 

DOUBLE  PRECISION  DSORT , UDS , VDS . UAS . VAS , RED . CD . 

. GS , RHOA , RHOD . RDS . NUS , HF 
C 

INTEGER  CDS 

C 

COMMON  / EQNMN/GS , RHOA , RHOD .  RDS  , NUS . HF 
C 

C  IN  UDS  = 

C  IN  VDS=DROPLET  VELOCITY  COMPONENTS. 

C  IN  UAS= 

C  IN  VAS -A  I R  VELOCITY  COMPONENTS. 

C  IN  CDS=PARAMETER  TO  DETERMINE  DRAG  COEFFICIENT  FORMULATION. 

C  OUT  RED  =  RELAT I VE  MOTION  REYNOLDS  NO. 

C  OUT  CD=DRAG  COEFFICIENT. 

C 

RED  =  DSORT ( ( UDS -UAS ) *  *  2+ ( VDS* VAS ) * *2 ) *  2 . DO*RDS/NUS 
I F ( CDS .  EQ .  2 )GOTO  300 

I F ( CDS . EQ . 1 . AND. RED. LE. 5. DO) GOTO  100 

C 

C  STEADY  STATE  DRAG  COEFFICIENT  OF  DROPLET  FOR  RED  <  5000 
C  ABRAHAM  (1970) 

CD=0. 2924D0M 1 +9 . 06D0/DSQRT ( RED ) )**2 
RETURN 

100  IF(RED.GE. 1 .D-2)G0T0  200 
C 

C  STEADY  STATE  STOKE ' S  DRAG  FOR  RED  <  0.01 
CD=24 .DO/RED 
RETURN 

C 

C  STEADY  STATE  DRAG  COEFFICIENT  FOR  0.01  <  RED  <  5  -  SARTOR 
C  AND  ABBOTT  ( 1975 ) 

200  CD  =  24 . DO/R  ED  +  2 . 2 DO 
RETURN 

C 

C  STEADY  STATE  DRAG  COEFFICIENT  -  LANGMUIR  &  BLODGETT  (1945) 

300  CD  =  24 . DO/RED+4  7  3DO/RED* +0 . 37DO+6 . 24D-3 *RED*  *0 . 38D0 

RETURN 
END 

C 

c 

SUBROUTINE  HIST(T.G) 

C 

C  WRITTEN  BY:  M.  OLESKIW  0N:8O1216  LAST  MOD  I F I  ED : 80 1 222 

C 

C  DETERMINES  VALUE  OF  INTEGRAL  IN  HISTORY  TERM  FOR  U  COMPONENT  EON 
C  REF:  BURDEN.  R.L..  U.D  FA  I  RES ,  &  A.C  REYNOLDS  (1978) 

C  NUMERICAL  ANALYSIS  P.  90  OA  297. B84 

C 
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DOUBLE  PRECISION  TAU3.TAU2.T AU 1 , T AUO ,  P  1 1 , P 10 , P2 1 ,P22, P20 , 
P33.P32.P31 . P30 , TO , T 1 . T2 . T3 . TS( 500 ) .F0.F1 , F 2 . F3 . DSQRT , DTS ( 6 ) . 

.  HT(2,6),T,A,B,C.D,F, AN( 2,6) ,P(2,745) .Z2.Z33.Z32.Z31 , Z30 , AA , BB 
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C 

INTEGER  U. I .  L.FF , E , MOD , JI .  J  J .  G , I , I M4 , I M3 , I M2 , I M 1 , 10. IP  1 
C 

COMMON  /LOC/TS.DTS. I , IM4 , I M3 , I M2 . I  Ml , 10, I P 1 / INTEG/AN , HT 

C 

C  IN  T=T IME  AT  PRESENT  TIME  STEP. 

C  IN  G=0: EXTRAPOLATE  HISTORY  TERM  SEQUENCE. 

C  IN  1 : CALCULATE  NEW  HISTORY  TERM. 

C 

C  STATEMENT  FNS .  TO  EVALUATE  PORTIONS  OF  THE  INTEGRAL. 

T AU3 ( A , B )  =  ( ( 5 . DO* A  *  *  3  +  6 . DO* A  * A  *T *8 . DO* A*  T  *  T+ 16 . DO*T*  *3 ) 

. * DSQRT ( T - A ) - ( 5 . DO  +  B *  * 3+6 . DO*B  +  B  * T  +  8  . DO*B* T  *  T  +  16 . DO*T  *  *3 ) 

. * DSQRT ( T-B ) ) *2 .DO/35  DO 

TAU2( A , B )  =  ( ( 3 . DO*  A  * A  +  4 . DO*A  *T  +  8 . DO*T  *T ) * DSQRT ( T - A ) 

. -( 3  D0  +  B*B  +  4 . DO*B  *  T  +  8 . DO*  T  *T ) * DSQRT (T-B ) )+2  .DO/ 15  DO 
T  AU 1 ( A . B )  =  ( ( 2 . DO*T+A ) +DSQRT (T-A)-(2. DO*T+B ) * DSQRT ( T-B ) ) *2 .00/3. DO 
TAUO( A,B)=2 . DO* (DSQRT (T-A) -DSQRT (T-B)  ) 

C 

C  STATEMENT  FNS.  TO  FIND  THE  TERMS  OF  THE  LAGRANGE  POLY.  FIT. 

PI 1(T0)=(F 1 -FO ) / ( T 1 -TO ) 

P10( T0)=( F0*T1-F 1*T0)/(T1-T0) 

22(A,B.C,F)=F/(A-B)/(A-C) 

P22(TO)=Z2(TO.T1 . T2 . FO )+Z2 ( T 1 .TO.T 2. FI )  +  Z2 ( T 2 , TO , T 1 ,F2) 

P2 1 ( T0)=-(T 1+T2 )* Z2( TO.T 1 . T2 . F0)-(T0+T2 ) • Z2 ( T 1 . TO , T 2 . F 1 ) 

. -( TO+T 1 )*Z2(T2.TO,T 1 ,F2) 

P20( TO)=T 1*T2*Z2(TO,T 1 ,T2,FO)+TO*T2*Z2(T 1 . TO. T2 . F 1 ) 

. +T0*  T 1*Z2(T2,TO.T 1  ,  F2  ) 

Z33(A,B.C.D,F)=F/(A-B)/(A~C)/(A-0) 

P33( TO )  =  Z33 ( TO. T 1.T2.T3.FO)+Z33(T1.TO.T2.T3.F1) 

. +  Z33( T2. TO. T 1 . T3 , F2 )  +  Z33( T3 . TO. T 1  . T2 . F3  ) 
Z32(A,B.C.D.F)=-(B+C+D)*F/(A-B)/(A-C)/(A-D) 

P32( TO)=Z32( TO.T 1 ,T2,T3.F0)+Z32(T 1 . TO . T 2 . T3 . F 1 ) 

. +  Z32 ( T  2 , TO. T 1 , T3, F2 )+Z32( T3.TO.T1.T2.F3) 

Z31 (A.B.C.D.F )=(B*C+B*D+C*D)*F/(A-B)/(A-C)/(A-D) 

P3 1 (T0)=Z3 1 ( TO. T1 . T2 . T3 . FO ) +  Z3 1 ( T 1 . TO . T 2 , T3 . F 1 ) 

. +Z3M  T2. TO.T 1 .T3.F2 ) +  Z3 1 <  T3 . TO, T 1 ,T2,F3) 

Z30( A.B.C.D.F )=-B*C*D*F/(A-B)/(A-C)/(A-D) 

P30( TO) =Z30( TO. T 1 , T2 , T3 , FO)+Z30( T1.TO.T2.T3.F1) 

. +  Z30( T2 . TO, T 1 , T3 , F2 )+Z30( T3.TO.T1.T2.F3) 

C 

I F ( G . EQ .  1  )GOTO  200 

C  EXTRAPOLATION  OF  HISTORY  TERM  SEQUENCE 
GOTO ( 140, 120. 100) . I 
TOrTS ( 1-3) 

T  1  =TS( 1-2) 

T  2 - T  S ( I-  1  ) 

T  3  =  T  S  (  I  )  . 

DO  110  J=1,2 
FO-HT ( J , I  M3 ) 

F 1 =HT ( J ,  I  M2 ) 

F  2  =HT ( J . IM 1  ) 

F  3  =  HT ( J , 10) 

HT ( J . IP  1 )=P33(T0)*T**3  +  P32<  TO ) *T * T  +  P3 1 ( TO ) *T  +  P30( TO ) 

110  CONTINUE 

RETURN 

C 

lOO  TO  =  TS(  1  ) 

T  1  =  TS( 2  ) 

T2=TS(3) 

DO  130  J=1,2 
FO=HT ( J , I M2 ) 

F 1 -HT ( J , I M 1 ) 

F2=HT( J. 10) 

HT( J. IP  1  )=P22< TO ) *  T  *  T  +  P2 1 ( TO) *T+P20( TO) 

130  CONTINUE 

RETURN 
C 

120  T0  =  T  S (  1  ) 

T 1 =T  S ( 2  ) 

DO  150  1 . 2 

FO=HT ( J , I M 1  ) 
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F 1 =HT ( J, 10 ) 

HT  ( J , IP  1 )  =P  1 1(T0)*T+P10(T0) 

150  CONTINUE 

RETURN 
C 

140  HT (  1 , IP  1 ) =0 . 00 

HT ( 2 , I P 1 ) =0 . DO 
RETURN 
C 

200  L=( I -4  )/2 *3+ 1 

HT ( 1 , I P 1 ) =0 . DO 
HT ( 2 , I P 1 ) =0 . DO 
GOTO ( 400 . 500 , 600 . 700 ) , I 
FF  =MOD (1.2) 

E  =  I -5+FF 

C  EVALUATE  INTEGRAL  UP  TO  LAST  SEVERAL  INTERVALS 
DO  210  J=1  .E.2 
AA  =  TS( J ) 

BB  =  TS ( J  +  2 ) 

JI=( J-1 )/2*3+1 
DO  220  J J= 1 . 2 

HT ( JU . I P 1 ) =HT ( JJ . I P 1 )+P( JJ. JI ) +  TAU2 ( AA , BB ) 

+  P( JJ, JI-M  )*TAU1(AA.BB)+P( JJ, JI +2 ) *T AUO( AA , BB ) 

220  CONTINUE 

210  CONTINUE 

IF( FF . EQ .  1  ) GOTO  600 

C  EVALUATE  INTEGRAL  FOR  LAST  4  INTERVALS 
C  USING  TWO  INTERVAL  PAIRS  (FOR  I  EVEN) 

700  T0=TS(I-3) 

T1=TS(I-2) 

T2=TS( I- 1 ) 

DO  710  J=1 ,2 
FO=AN( J, IM3 ) 

F 1 =AN( J , I  M2  ) 

F2=AN( J , I M 1 ) 

C  FIT  A  2ND  ORDER  LAGRANGE  POLYNOMIAL 
P( J.L)=P22(TO) 

P( J.L+1 ) =P2 1 ( TO ) 

P( J , L+2 ) =P20( TO ) 

HT  ( J , I P 1 ) =HT ( J , IP1 )+P( J.L)*TAU2(TO.T2)+P( J.L+1 ) *TAU1 ( TO, T2 ) 
+P< J.L+2)*TAUO(TO,T2) 

710  CONTINUE 

C  FOR  THE  SECOND  PAIR  OF  THE  SET 
C  (OR  FOR  THE  VERY  FIRST  PAIR  OF  INTERVALS) 

500  T0=TS(I"1) 

T 1 =TS( I ) 

T2-TS(I+1 ) 

DO  720  J=1 ,2 
F0=AN( J , IM1 ) 

F  1  = AN( J . 10) 

F2=AN( J , IP1 ) 

HT ( J  f I P 1 ) *HT ( J , I P 1 )+P22(TO)*TAU2(TO,T2)  +  P21(TO)*TAU1(TO,T2) 
+  P20(T0)*T  AUO( TO , V2) 

720  CONTINUE 

RETURN 
C 

C  EVALUATE  INTEGRAL  FOR  LAST  3  INTERVALS  (FOR  I  ODD) 

600  T0=  TS ( I ~  2 ) 

T1  =  TS(I'1  ) 

T2  =  TS( I  ) 

T3=T$( 1+1 ) 

DO  610  J=1 .2 
F0= AN( J , I M2 ) 

F 1 = AN( J , I M 1 ) 

F 2  - AN( J ,  10) 

F  3  =  AN( J , IP1 ) 

HT( J. IP1 ) -HT ( U , IP1 )+P33(TO)*TAU3(TO,T3)+P32(TO)*TAU2(TO.T3) 
+P31(TO)*TAU1(TO.T3 )+P30( TO) #TAUO( TO , T3 ) 

610  CONTINUE 

RETURN 
C 

C  EVALUATE  INTEGRAL  FOR  THE  FIRST  INTERVAL 
400  TO=TS( 1 ) 

T  1  —  T  S  (  2  ) 
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DO  410  d  = 1 . 2 
F0=AN( J. 10) 

F 1 =  AN( d , I P 1  ) 

HT ( J . I P 1 ) =HT ( J , I P 1 )+P1 1(T0)*TAU1(T0,T1 ) +P 10( TO ) *TAUO( TO . T 1 ) 
410  CONTINUE 

RETURN 
END 

C 

C 

SUBROUTINE  RKF4 ( EON . CDS . EPS ) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:  800227  LAST  MOD I F I  ED : 80 1 227 

C 

C  INTEGRATE  THE  DROPLET  EONS.  OF  MOTION  (IN  X  AND  Y)  USING 
C  THE  4TH  ORDER  RUNGE -KUTTA-FEHLBERG  TECHNIQUE. 

C  REF:  BURDEN,  R  .  L  .  ,  d.D.  FAIRES,  &  A.C.  REYNOLDS  (1978), 

C  NUMERICAL  ANALYSIS,  P . 254 ,  QA  297 . B84 

C 

DOUBLE  PRECISION  EPS , XDS ( 6 ) , UDS ( 6 ) , AN( 2 , 6 ) . YDS ( 6 ) , 

. VDS(6) , HT (2,6) . DTS ( 6 ) ,UAS(6) , VAS(6) , RED(6 ) ,CD,RE, 

. K 1 ,K2 . K3 , K4 , K5 . K6 . L 1 , L2 . L3 , L4 , L5 , L6 .Ml . M2 . M3 , M4 , M5 , M6 , 

. N 1 ,N2 ,N3,N4 .N5.N6.UA, VA ,RMAX ,DMAX1 .DMIN.0MIN1 , XDEL , YDEL , 

UDEL. VDEL.DABS.XR. YR.UR.VR.XT, YT.UT. VT, 

. XD , YD tUD , VD , CC 1 , CC2 , C3 , C4 , C5 , C6 .C7.C8.C9, CIO, Cl  1 .C12.C13, 
.C14.C15.C16.C17 . C 18 ,C 19, C20. C2 1 , C22 , C23 , C24 , TS( 500) 

DOUBLE  PRECISION  DMINP 

C 

INTEGER  EON, CDS, I  , I M4 , I M3 , I  M2 , I M 1 , 1 0 . I P 1 . K 

C 

COMMON  /PV/XDS , YDS , UDS , VOS/ I NTEG/ AN , HT 
. /LOC/TS.DTS, I , IM4 , IM3, I M2, IM1 , 10, IP  1 
/REL/UAS, VAS.RED.CD 

. / RKFM/CC 1 ,CC2 .C3.C4 .C5.C6.C7 , C8 , C9 , C 10 . C 1 1 ,C12,C13,C14, 
.C15,C16,C17,C18,C19.C20,C21 .C22.C23.C24 
C 

C  IN  EON=DENOTES  PORTION  OF  TOTAL  SYSTEM  OF  EQUATIONS  TO  BE  SOLVED. 
C  IN  CDS=T YPE  OF  DRAG  COEFFICIENT  TO  BE  USED. 

C  IN  E PS  =  LOC A L  ERROR  PARAMETER. 

C 

IF (DABS (DMINP ) . L  T .  1  . D-70  )DMI NP= 1.01  DO 
100  TS (  1  +  1  )  =  TS ( I  )+DTS( IO) 

K 1 =DT  S ( 10 ) *UDS (  10) 

L 1 =DTS ( 10 ) * VDS ( 10) 

M 1 =DTS ( 10 ) *  AN(  1.10) 

N 1 =DTS ( 1 0  )  +  AN ( 2,10) 

XD=XDS ( I0)+CC1 *  K 1 
YD  =  YDS (  IO)*CC  1  * L  1 
UD=UDS( IO)+CC 1 "Ml 
VD - VDS ( IO)+CC 1  *N1 
CALL  AIRVEL(XD. YD.UA.VA.4) 

CALL  DRAG(UD. VD , UA . V A , CDS , RE . CD ) 

C 

K2=DIS( 10) "UD 
L2^DTS( 10 ) * VD 

CALL  ACCN(UO . VO . UA , V A , RE , CD , EON, TS( I )+DTS( I0)/4 . DO.O) 

M2=DTS( 10) "AN( 1 , IP1 ) 

N2=DTS(  I0)'*AN(2.  I  P  1  ) 

XD  =  XDSI  I0)+CC2"K 1+C3*K2 
YD  =  YDS  (  I0)+CC2*L  1  +-C3^L2 
UD=UD5( IO)+CC2"M1+C3"M2 
VD  =  VDS (  IO)+CC2*N1+C3*N2 
CALL  AIRVEL(XD. YD.UA.VA.4) 

CALL  DRAG(UD, VD .UA , VA .CDS . RE .CD ) 

C 

K3=DTS( 10 ) *UD 
L3=DTS( 10 ) * VD 

CALL  ACCN( UD , VD , UA , VA , RE . CD, EON. TS< I )+DTS( 10) *3 . 75D- 1,0) 

M3  =  DT  S (  10 ) *  AN(  1  , IP  1  ) 

N3=DTS( 10) *AN( 2 . IP1 ) 

XD  =  XDS (  I0)+C4*K 1 -C5*K2  +  C6*K3 
YD=YDS( I0)+C4*L 1 ~C5*L2+C6*L3 
UD=UDS ( I0)+C4 "Ml -C5"M2+C6"M3 
VD - VDS ( I0)  +  C4 *  N 1  -C5"N2  +  C6"N3 
CALL  AIRVEL(XD, YD.UA.VA.4  ) 
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CALL  DRAG(  LID  ,  VD  ,  UA  ,  VA  ,  CDS  ,  RE  ,  CD  ) 

C 

K4=DTS ( 10 ) *UD 
L4=DTS( 10 ) * VD 

CALL  ACCN(UD, VD.UA, VA.RE .CD, EQN.TS( I )+ 1 2 . DO/ 1 3 . DO 
. *DTS  ( 10) ,0) 

M4=DTS( I0)*AN( 1 , IP1 ) 

N4  =  DTS ( 10 ) * AN( 2 , I P 1 ) 

XD=XDS( I0)+C7*K 1-C8*K2+C9*K3-C10*K4 
YD  =  YDS ( I0)+C7*L1-C8*L2+C9*L3-C10*L4 
UD=UDS( I0)+C7+M1 -C8*M2+C9*M3-C10*M4 
VD= VDS ( I0)+C7*N1-C8*N2+C9*N3-C10+N4 
CALL  AIRVEL(XD, YD , UA , VA ,4 ) 

CALL  DRAG( UD , VD.UA, VA.CDS, RE, CD ) 

C 

K5=DTS( 10 ) *UD 
L5=DTS ( 10 ) * VD 

CALL  ACCN(UD, VD.UA . VA , RE . CD . EON . TS ( I + 1 ) .0) 

M5  =  DT  S ( 10 ) *  AN(  1 , IP1  ) 

N5=DTS( I0)*AN(2, IP1 ) 

XD  =  XDS( 10) "Cl  1 *K1+C12*K2-C13*K3+C14*K4-C15*K5 
YD- YDS( 10) -Cl  1 +L1+C12*L2'C13*L3+C14*L4-C15*L5 
UD  =  UDS( 10) -Cl  1 +M1+C12*M2~C13*M3+C14*M4-C15*M5 
VD  =  VDS  (  1 0  )  -C  1  1  *N1+C12*N2-C13*N3+C14*N4-C15*N5 
CA!  L  AIRVEL(XD. YD.UA.VA. 4) 

CALL  DRAG(UD, VD.UA, VA , CDS . RE , CD ) 

C 

K6=DTS( I0)*U0 
L6=DTS( 10 ) * VD 

CALL  ACCN(UD. VD.UA. VA.RE .CD. EQN,TS( I )+DTS ( 10 )/2 . DO . 0 ) 

M6=DTS ( 10 ) *  AN (  1  . IP1  ) 

N6  =  DTS( 10)*  AN ( 2 , IP1 ) 

C 

C  NEW  DROPLET  POSITION  AT  1+1 

XDS ( IP1 )  =  XDS  ( I0)+C16*K1+C17*K3+C18*K4-C19*K5 
YDS ( I P  1  )  =  YDS ( I O  )  +C 16*L1+C17*L3+C18*L4-C19*L5 
C  NEW  DROPLET  VELOCITY  AT  1+1 

UDS ( IP1  )  =  UDS  (  IO)+C 16*M1+C17*M3+C 1 8*M4-C 19*M5 
VDS( IP  1 ) = VDS ( I0)+C16*N1+C17*N3+C18*N4-C19*N5 
C 

C  5TH  ORDER  ESTIMATE  OF  POSITION  AND  VELOCITY 

XT  =  XD$ (  10 ) +C20*K 1 +C2 1 *K3+C22 *K4 -C23 * K5  +  C24 *K6 
YT  =  YDS ( 10 )+C20*L 1+C2 1  * L3  +  C2 2 ♦ L4 -C23 ♦ L5  +  C24 * L6 
U  T - UDS ( 10 ) +  C20  *M 1 +C2 1  * M3 +C22 ♦ M4 -C23  *  M5  +  C24 * M6 
VT  =  VDS  (  I0)+C20*N1+C2 1  *N3  +  C22  *N4 -C23  *N5*C24  *N6 
C 

C  DETERMINE  DEFFERENCES  IN  4TH  AND  5TH  ORDER  ESTIMATES. 

XR  =  ( XT- XDS ( IP1  )  )/DTS(  10) 

I F ( DABS ( XR  )  . L  T .  1 . 0~ 70 )XR  = t . D- 70 
YR  =  ( YT-YDS(  IP  1  )  )/DTS( 10) 

I F ( DABS ( YR  )  .  LT  1 . D - 70 ) YR = 1 . D - 70 
UR  = ( UT -UDS (  IP1  )  )/DTS(  10) 

IF (DABS (UR)  LT . 1  D~70)UR= 1 .0-70 
VR “ ( V  T - VDS ( IP1 ) )/DTS( 10) 

I F ( DABS ( VR  )  . LT .  1  D-70)VR= 1  .D-70 
C  CALCULATE  STEP  SIZE  ADUST  I NG  FACTORS. 

XDEL  =  ( EPS/DABS ( XR  )  )** . 25DO 
YDE L  =  ( E PS/DABS ( YR  )  )*♦  . 25DO 
UDEL  =  ( EPS/DABS (UR  )  )* ♦  25DO 
VDEL  =  ( EPS/DABS ( VR)  )* ♦  25DO 
C  ADJUST  FOR  LEAST  PRECISE  EON. 

DMIN  =  DMI N 1 ( XDEL . YDEL . UDE  L . VDE  L ) 

RMAX  =  DMAX 1 ( DABS ( XR ) . DABS( YR) .DABS (UR)  , DABS ( VR  )  ) 

K  =  I O 

I F ( RMAX . LF .  EPS  )K  =  I P 1 
I F ( DMI NP . L  T  1  . DO  )GOTO  200 

I F ( DMI N . LT . 1  D0)DTS(K)=0.9D0*DMIN+DTS( 10) 

I F ( DMI N . GE . 1 1  D0)DTS(K)=1  800*DTS( 10) 

IF (DM IN. GE  1  DO  AND  DMI N . LT . 1 1  DO )DTS< K ) = ( (DMIN- 1 . DO)/ 10  D0+ 1 . DO ) 
*0  T  S (  r0)*O.9DO 
GOTO  210 

200  IF (DMIN. LE  O  5D0 )DTS < K ) =0 . 5D0*DTS ( 10) 

IF (DMIN  GT  1  D0)DTS(K)=( (DMIN- 1  DO)/ 10.D0+1 . DO ) * 0 . 9D0*DT S ( 10) 

I F ( DMIN . GT  0.5D0  AND. DMIN. LE .  1  . DO )DTS( K ) =DMI N*0 . 9D0*DT S ( 10) 


3018 

3019 

3020 
302  1 
3022 

302  3 

3024 

3025 

3026 

3027 

3028 

3029 

3030 

303  1 

3032 

3033 

3034 

3035 

3036 

3037 

3038 

3039 

3040 

304  1 
304  2 

3043 

3044 

3045 

3046 

304  7 

3048 

3049 

3050 

3051 

305  2 

3053 

3054 

3055 

3056 

3057 

3058 

3059 

3060 

306  1 

3062 

3063 

3064 

3065 

3066 

3067 

3068 

3069 

3070 

307  1 

3072 

3073 

3074 

3075 

3076 

3077 

3078 

3079 

3080 

308  1 
308  2 

3083 

3084 

3085 

3086 

308  7 

3088 

3089 

3090 

309  1 


210  DMI NP=DMIN 

I  F  ( RMAX  GT . EPS  )G0T0  100 

C 

C  NEW  ACCELERATIONS  AT  1+1 

CALL  AIRVEL( XDS( IP  1 ) #Y0S(IP1 )  ,UAS( IP  1  )  . VAS( IP  1 ) , 5 ) 

CALL  DRAG( UDS  l IP  1 ) . VDS( IP1 ) ,UAS( IP  1 ) , VAS( IP1 ) , CDS , RED( IP1 ) .CD) 
CALL  ACCN ( UDS ( IP  1 ) . VDS( I P 1 ) ,UAS( IP  1  )  , VAS( IP1 ) . RED( I P 1 ) . 

CD. EON. TS( 1+1 ).0) 

IF( EON.NE . 2  )RE TURN 

CALL  ACCN( UOS ( IP1 ) ,  VOS( IP  1 ) ,UAS( IP1 ) ,VAS( IP1 ) ,REO( IP1 ) , 

CD , EQN . TS ( 1  +  1  ),  1  ) 

RETURN 

END 

C 

C 

SUBROUTINE  RK4(E0N,CDS) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:  790926  LAST  MOD I F I  ED : 80 1 223 
C 

C  INTEGRATE  THE  DROPLET  EONS  OF  MOTION  (IN  X  AND  Y)  USING  THE  4TH 
C  ORDER  RUNGE-KUTTA  TECHNIQUE. 

C  REF:  BURDEN. R.L.,  U.D.  FA  I  RES .  8  A.C.  REYNOLDS  (  1978).  NUMERICAL 
C  ANALYSIS  P.  281  QA  297. B84 

C 

DOUBLE  PRECISION  K 1 . L 1 . K2 . L2 , K3 . L3 . K4 . L4 . DTS ( 6 ) , XDS ( 6 ) . UDS( 6 ) . 
. YDS (6 ) . VDS ( 6 ) , AN( 2,6 ) ,HT(2 .6) . 

Ml , M2 , M3 . M4 . N 1 , N2 .N3.N4.U1 .U2.U3.V1 .V2,V3,CD,RE.RED(6), 
VAS(6),UAS(6).TS(500) 

C 

INTEGER  I , EON. I M4 . I M3 . I M2 , I M 1 , I O . I P 1 .CDS 
C 

COMMON  / I  NT  EG/ AN . HT/PV/XDS . YDS . UDS . VDS 
. /LOC/TS .DTS,  I , IM4 , I M3. I M2 . IM1 . 10, IP  1 
. /REL/UAS, VAS.RED.CD 
C 

C  IN  DTS=NON -DIMENSIONAL  TIME  STEP 
C  IN  1=  PRESENT  INDEX  OF  VECTORS  XDS. UDS.... 

C  IN  EQN=  CHOICE  OF  TERMS  USED  IN  RHS  OF  ODE 
C 

TS ( 1+1  )  - TS ( I  ) +DTS ( 10) 

K  1  =D  TS ( 10) *UDS ( 10) 

L  1  =DTS ( 10) *VDS( 10) 

M 1  -DTS ( 10) *AN(  1 . 10) 

N  1  =0  T  S (  10)*  AN ( 2 , 10) 

CALL  A  I RVE  L ( XDS ( I O ) +K 1 /2 . DO . YDS ( IO)  +  L 1/2 . DO.U 1 . VI . 4 ) 

CALL  DRAG( UDS ( I0)+M1/2 . DO. VDS ( 10 )+N 1/2 . DO. U 1 . V 1 , CDS . RE , CD ) 

C 

K2  -DT  S ( 1 0 ) * ( UDS (  I0)+M1/2.D0) 

L2  =DTS ( 10) M  VDS (  I0)+N1/2 . DO) 

CALL  ACCN ( UOS ( 10 ) +M 1 /2 . DO. VDS( I O ) +N 1 / 2 . DO . U 1 . V 1 . RE .CD. EON, 
TSID.O) 

M2  =DTS ( 10) *AN(  1 , IP1  ) 

N2  =  DTS ( 10)*  AN( 2 , I P 1  ) 

CALL  A  I RVE L ( XDS ( I O ) +K2/2 • DO , YDS ( IO)  +  L2/2  D0.U2 . V2 ,4 ) 

CALL  DRAG( UDS ( I 0 ) +M 1 /2 . DO , VDS ( 10 )+N 1 /2 . DO . U2 . V2 , CDS . RE . CD ) 

C 

K3  =  DT S ( 10) *(UDS( IO)+M2/2  00) 

L  3  =  DT  S ( I O ) * ( VDS (  I O ) +N2/2 . DO ) 

CALL  ACCN ( UDS ( I O ) +M2/2 . DO . VDS ( 10 ) +N2/2 . DO . U2 , V2 . RE . CD , EON , 

. T S ( I  )  +DTS ( I0)/2 .00,0) 

M3 - DT  S (  IO)*AN(  1  ,  IP1  ) 

N3-DTS( 10 ) * AN( 2 , IP1 ) 

CALL  A  I RVE  L ( XDS (  I0)+K3. YDS( I O )  +  L3 , U3 , V3 , 4 ) 

CALL  DRAG( UDS ( I0)+M3, VDS( I 0 ) +N3 , U3 , V3 , CDS , RE .CD) 

C 

K4  =  D  T  S ( I O ) * ( UDS (  I0)+M3) 

L4=DTS( 10) *( VDS( I0)+N3) 

CALL  ACCN ( UDS ( I0)+M3.VDS( 10 ) +N3 , U3 , V3 , RE . CD . EON . 

T S ( I  )  +DTS ( I0)/2  DO. O) 

M4  =DT  S (  10  )  *  AN(  1  .  IP1  ) 

N4  =  DT  S( 10 ) *  AN( 2 .  IP1  ) 

C 

C  NEW  DROPLET  POSITION  AT  1+1 

XDS (  IP1  )  =  XDS (  IO)  +  (K1+2.DO*K2  +  2.DO*K3+K4)/6.DO 
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YDS ( I P 1 ) =  YDS ( I0)+(L 1+2  D0*L2+2.D0*L3+L4)/6 .DO 
C  NEW  VELOCITIES  AT  1+1 

UD$( IP  1  ) =UDS ( I0)  +  (M1  +  2 . DO* M2 +  2 . DO*M3+M4 )/6 . DO 
VDS ( IP  1 )  = VDS( I0)  +  ( N1  +  2 . DO*N2+2 . D0*N3+N4 ) /6 . DO 
C  NEW  ACCELERATIONS  AT  1+1 

CALL  A I RVE  L ( XDS  (  IP1 ) , YDS( IP1 )  tUAS( IP1 )  ,VAS( IP1  ) .5) 

CALL  DRAG (  LIDS (  IP1  )  .  VDS(  I  P  1  )  .  UAS<  IP  1  )  ,  VAS(  IP1  )  .  CDS  .  RED(  I P  1  )  ,CD) 
CALL  ACCN( UDS ( IP1 ) . VDS( IP1  )  .UAS( IP1  )  , VAS( IP1  )  , RED( IP1 ) .CD. EON. 
TS ( 1+1 ) .0) 

DTS< IP1  )=DTS(  10) 

I  F  ( EON  NE  .  2 )RETURN 
C 

CALL  ACCN ( UDS ( IP1  )  ,VDS( IP1  )  ,UAS( IP1  )  .VAS( IP1  ) , RED( I P 1 ) .CD. EON. 
TS( 1+ 1 ) . 1 ) 

RETURN 

END 

C 

C 

SUBROUTINE  PC4(EQN,CDS) 

C 

C  WRITTEN  BY:  M.  OLESKIW  ON:  800122  LAST  MODIFIED:  801223 

C 

C  INTEGRATE  EONS.  OF  MOTION  USING  THE  4TH  ORDER  PREDICTOR- 
C  CORRECTOR  METHOD 

C  REF:  BURDEN.  R.L.,  J.D.  F  A I  RES  t  &  A.C.  REYNOLDS  (  1978), 

C  NUMERICAL  ANALYSIS  OA  297. B84  P  .  266 

C  HAMMING.  R.W.  (1973).  NUMERICAL  METHODS  FOR  SCIENTISTS  & 

C  ENGINEERS.  2ND  ED..  OA  297. H28  CHAPS.  22  S  23 

C 

DOUBLE  PRECISION  XDS ( 6  )  . UDS ( 6 ) . AN( 2 , 6 ) , HT ( 2 , 6 ) , YDS ( 6 ) , 

. VDS (6 ) , AO, A  1 . A2.B0.B1 ,B2 ,B3. 

•  CO. C 1 . C2 ,DM1 . DO.D 1 .02 .UP  I , UCI , VP  I . VC  I , MUAS . MVAS, 

PUDS .DTS(6  )  , PVDS . MUOS , MVDS . CUDS , CVDS . UDSP 1 . VDSP 1 
, F  MU , FMV.UST.VST.ERI . ER2 . PXDS . P YDS . MXDS . MYDS . CXDS . C YDS 
, UAS ( 6 ) , VAS ( 6  )  . RED ( 6  )  , XP I , XCI  . YPI , YCI . RE . CD . TS ( 500 ) 

C 

INTEGER  I . EON,  IM4 .  I  M3.  I M2 . IM1 . 10. IP  1 .CDS 
C 

COMMON/ I  NT EG/ AN . HT/PV/XDS . YDS . UDS . VDS 

/PCM/ AO. A  1  . A 2 .BO, B 1 .B2 , B3 . CO. C 1 .C2 . DM1 . D0.D1 ,D2 . 

UP  I  .UCI  , VP  I . VC  I , ER 1 , ER2 , XPI , XCI , YPI , YCI ,UST , VST 
. /LOC/TS.DTS .  I . IM4,  I  M3. I  M2. IM1 . 10.  I P 1 
./REL/UAS.VAS.RED.CD 
C 

C  IN  EQN  =  CHOICE  OF  TERMS  USED  IN  RHS  OF  ODE 
C  IN  CDS=TYPE  OF  DRAG  COEFFICIENT  TO  BE  USED. 

C 

TS(  1  +  1  )  =  T S ( I  )  +  DT S ( 10) 

C 

C  THE  PREDICTOR 

PXDS  =  A0*XDS(I0)  +  A1 *XDS( IM1 )+A2*XDS(  I  M2 ) 

♦DTS(  10) MBO*UDS( I0)+B1 *UDS<  IM1  )+B2*UDS(  I  M2 )+B3*UDS(  I M3 ) ) 

P  YDS  =  AO  *  YDS (  10)  +  A  1  * YDS( IM1)+A2*YDS(IM2) 

♦DTS(  10) M BO* VDS ( 10 ) +B 1 *VDS( IM1 )+B2*VDS< 1  M2 ) +B3 * VDS ( I M3 ) ) 

PUDS -AO* UDS ( I O ) * A  1 *UDS( I M 1 )+A2*UDS<  I  M2 ) 

+  DTS(  10) ♦ (BO*AN(  1 ,  IO)+B 1  *AN(  1 .  IM1 )+B2*AN(  1  .  I  M2  )+B3*AN(  1 .  IM3 )  ) 
PVDS  =  AO* VDS (  I O )  +  A  1  * VDS( IM1  )+A2*VDS<  I  M2 ) 

*DTS(  I0)*(B0*AN(2.  I0)+B1*AN(2,  I M  1  )+B2*AN(2,  I  M2  )  +B3  *  AN(  2  ,  IM3)) 

C 

C  MODIFICATION  OF  THE  PREDICTOR 
MXDS  =  PXDS-ER1 * ( XPI-XCI  ) 

MYDS  =  PYDS-ER  IMYPI-'iCI  ) 

MUDS -PUDS  - ER 1 ♦ ( UP  I -UCI  ) 

MVDS "PVDS - ER  1  *( VP  I  VC  I  ) 

CALL  A  I RVE L( MXDS .MVDS. MUAS. MVAS. 4 ) 

CALL  DRAG ( MUDS , MVDS . MUAS , MVAS , CDS , RE , CD ) 

CALL  ACCN( MUDS . MVDS , MUAS , MVAS . RE , CD . EON. TS( I + 1 ) . 0 ) 

FMU-AN1  1  .  IP  1  ) 

F  MV  =  AN( 2 . IP1  ) 

C 

C  THE  CORRECTOR 

CXDS  =  CO  *  XDS (  IO)  +  C 1  *XDS( IM1  )+C2*XDS(  I M2 ) 

+  DTS<  10) M DM1  * MUDS  +DO*UDS<  I0)+D1 *UDS( I M 1 ) +D2 *UDS (  ] M2 ) ) 

C Y  DS -CO*  YDS (  I O ) +C 1  *  YDS (  IM1  )  +  C2*YDS(  I  M2) 


3166 

3167 

3168 

3169 
3  170 

3171 

3172 

3173 

3174 
3  1  75 
3176 
3  177 
3  178 
3179 
3  180 
3181 
3  182 
3183 
3  184 
3  185 
3  186 
3  187 
3188 
3  189 
3  190 
3191 
3  192 
3193 
3  194 
3  195 
3  196 
3  197 
3  198 
3  199 

3200 

3201 

3202 
END  OF 


. +DTS ( I0)*(DM1*MVDS+D0*VDS( I0)+D1*VDS( IM1 )+D2*VDS( I M2) ) 
CUDS=CO*UDS( I0)+C1*U0S( IM1 )+C2*U0S( IM2) 

.  +DTS(  I0)*(DM1*FMU+D0*AN(  1  ,  IO)+01*AN(  1  ,  IM1  )+D2*AN<  1  .  I  M2  )  ) 
CVDS-C0*VDS( I0)+C1*VDS( IM1 )+C2*VDS( IM2) 

.  +  DTS  (  I0)*(DM1+FMV+D0*AN(2,  I0)+D1*AN(2,  IM 1 )+D2*AN( 2 ,  I  M2)  ) 

C 

C  FINAL  VALUES 

XDS ( I P 1  )=CXDS+ER2*(PXDS-CXDS) 

YDS ( IP1  )=CYDS  +  ER2*(PYDS-CYDS) 

UDS  (  IP  1  )  =.CUDS  +  ER2  +  (  PUDS -CUDS  ) 

VDS( IP1 )=CVDS+ER2*( PVDS-CVDS ) 

C 

C  NEW  VALUES  FOR  ACCELERATION  AT  1+1 

CALL  A I RVEL ( XDS ( IP1 ) , YDS( I P 1 ) , UAS ( I P 1 ) , VAS ( I P 1 )  ,5  ) 

CALL  DRAG( UDS ( I P 1 ) . VDS( IP1 ) ,UAS( IP1 ) , VAS( IP1 ) .CDS, RED ( IP1 ) .CD) 
CALL  ACCN( UDS ( I P 1 ) . VDS( I P 1 ) , UAS ( I P 1 ) , VAS ( I P 1 ) , RED( IP  1 ) , CD , EQN, 
-TS(I-M).O) 

UDSP 1 =  AN(  1 .  IP1  ) 

VDSP 1 =  AN ( 2 , IP1  ) 

I F ( EQN . NE . 2 )GOTO  100 

CALL  ACCN ( UDS ( I P 1  ) , VDS ( I P 1 ) . UAS ( I P 1 ) . VAS ( I P 1 ) . RED ( I P 1 ) . CD , EON . 
. T  S ( 1  +  1 ) .  1) 

C 

C  CALCULATE  STABILITY  INDICES 

100  UST= ( F MU -UDSP 1 ) / ( MUDS -UDS ( I P 1 ) ) 

VST  =  (FMV- VDSP 1  )/ ( MVDS- VDS( IP  1 ) ) 

XPI =PXDS 
XCI =CXDS 
YPI =P YDS 
YC I =C YDS 
UPI=PUOS 
UCI =CUDS 
VPI=PVDS 
VC  I -C VDS 

DTS ( IP1 ) =DT  S ( 10) 

RETURN 

END 

FILE 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC 
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